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TURBULENT CONVECTION IN STELLAR INTERIORS. I. HYDRODYNAMIC SIMULATION 
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ABSTRACT 

We describe the results of three-dimensional (3D) numerical simulations designed to study turbu- 
lent convection in the stellar interiors, and compare them to stellar mixing-length theory (MLT). 
Simulations in 2D are significantly difi^erent from 3D, both in terms of flow morphology and velocity 
amplitude. Convective mixing regions are better predicted using a dynamic boundary condition based 
on the bulk Richardson number than by purely local, static criteria like Schwarzschild or Ledoux. 
MLT gives a good description of the velocity scale and temperature gradient for a mixing length 
of ~ l.lHp for shell convection, however there are other important effects that it does not capture, 
mostly related to the dynamical motion of the boundaries between convective and nonconvective re- 
gions. There is asymmetry between up and down flows, so the net kinetic energy flux is not zero. The 
motion of convective boundaries is a source of gravity waves; this is a necessary consequence of the 
deceleration of convective plumes. Convective "overshooting" is best described as an elastic response 
by the convective boundary, rather than ballistic penetration of the stable layers by turbulent eddies. 
The convective boundaries are rife with internal and interfacial wave motions, and a variety of insta- 
bilities arise which induce mixing through in process best described as turbulent entrainment. We find 
that the rate at which material entrainment proceeds at the boundaries is consistent with analogous 
laboratory experiments as well as simulation and observation of terrestrial atmospheric mixing. In 
particular, the normalized entrainment rate E—ue/cth, is well described by a power law dependance 
on the bulk Richardson number Ris = AbL/ajj for the conditions studied, 20 < Ris ^ 420. We 
find E = ARi]^", with best flt values, log A = 0.027 ± 0.38, and n = 1.05 ± 0.21. We discuss the 
applicability of these results to stellar evolution calculations. 

Subject headings: stars: evolution - stars: nucleosynthesis - massive stars - hydrodynamics - convection 
- g- modes 
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1. INTRODUCTION 

We have simulated three-dimensional (3D), turbulent, 
thermally-relaxed, nearly adiabatic convection (high 
Peclet number). Such flow is relevant to deep convec- 
tive regions in stars (i.e., to most stellar mass which is 
convective, but not mildly sub-photospheric and surface 
regions). We simulate oxygen shell burning on its natu- 
ral time scale, and core hydrogen burning driven at 10 
times its natural rate. The simulations develop a robust 
quasi-steady behavior in a statistical sense, with signif- 
icant intermittency. We analyze this statistical behav- 
ior quantitatively, and compare it to predictions of as- 
trophysical mixing length theory (Bohm-Vitense 1958). 
Mixing-length Theory (MLT) gives a good representa- 
tion of many aspects of convection, but omits others (es- 
pecially wave generation and mass entrainment) which 
are related to the dynamical behavior of stably stratified 
layers adjacent to the convection. 

In Section 2 we briefly summarize some results of the 
study of turbulent entrainment in geophysics, to prepare 
the reader for its appearance in our astrophysical sim- 
ulations. This process is not included in the standard 
approach to stellar evolution (Cox & Guili 1968; Clayton 
1983; Kippenhahn & Weigert 1990; Hansen & Kawaler 
1994). In Section 3 we discuss our numerical and theo- 
retical tools. In Section 4 we present our simulations of 
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oxygen shell burning, which attain a thermal steady state 
(this is possible because of the rapidity of nuclear heat- 
ing and neutrino cooling). In Section 5 we discuss a less 
advanced burning stage, core hydrogen burning, which 
we are able to examine with the use of an artiflcially en- 
hanced hydrogen burning rate (by a factor of ten). We 
find that the behavior is similar to the oxygen burning 
shell, suggesting that our results may have broad appli- 
cation for stellar evolution. In Section 6 we compare our 
results to the assumptions of MLT, and in Section 7 we 
show that our results lead to a simple model of turbulent 
entrainment, an effect not in MLT nor in standard stellar 
evolutionary calculations. 

This paper is the first in a series. In subsequent pa- 
pers, we incorporate the "empirical" convection model 
developed in this paper into the TYCHO stellar evolu- 
tion code (Young & Arnett 2005) and begin to assess its 
influence on stellar evolution, on nucleosynthetic yields, 
and on the structure of supernova progenitors. 

2. TURBULENT ENTRAINMENT 

The presence of a turbulent layer contiguous with a sta- 
bly stratified layer is common in both astrophysical and 
geophysical flows. Turbulence in a stratified media is of- 
ten sustained by strong shear flows or thermal convection 
and bound by a stabilizing density interface. Over time, 
the turbulent layer "diffuses" into the stable layer and 
the density interface recedes, thus increasing the size of 
the mixed, turbulent region. The basic features of this 
turbulent entrainment problem are illustrated in Figure 
1. The rate at which the density interface recedes into 
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the stable layer ue = dri/dt is called the entrainment 
rate, and its depcndance on the parameters character- 
izing the turbulent and the stable layers has been the 
subject of numerous experimental and theoretical stud- 
ies. It is generally ignored in stellar evolutionary studies. 

Experimental studies have mostly been of "mixing 
box" type which involves a tank of fluid with a turbu- 
lent layer and a density stratified layer. The turbulence 
is generated by thermal convection or an oscillating wire 
mesh, and density stratification imposed by either a so- 
lute or thermal gradient (Turner 1980). Complementary 
to these shear-free mixing box models arc shear drivcn- 
models. Shear-driven turbulence experiments involve ei- 
ther a recirculation track which propels one layer of fluid 
above a stationary layer, or a rotating plate in contact 
with the fluid that drives a circulation in the upper layer. 
Shear instabilities sustain a turbulent mixed layer in the 
overlying fluid which then entrains fluid from the lower, 
stationary layer (Kantha et al. 1977; Strang & Fernando 
2001). In all of these laboratory experiments, a variety of 
flow visualization techniques are used to study both the 
overall entrainment rate ue and the physical mechanisms 
which underly the entrainment process. 

One of the primary conclusions of these studies is that 
the entrainment rate depends on a Richardson number, 
which is a dimensionless measure of the "stiff'ness" of 
the boundary relative to the strength of the turbulence. 
In shear-free turbulent entrainment the bulk Richardson 
number. 



is most commonly studied. Here, A6 is the buoyancy 
jump across the interface, cr is the r.m.s. turbulence ve- 
locity adjacent the interface, and i is a length scale for 
the turbulent motions often taken to be the horizontal 
integral scale of the turbulence at the interface. The rel- 
ative buoyancy is deflned by the integral, 

r 

b{r) = j N^dr (2) 
where N is the buoyancy frequency defined by. 

The entrainment coefficient E is the interface migration 
speed Ue normalized by the R.M.S. tur blent velocity at 

the interface E = ue/(t, and is generally found to obey 
a power law dependanco on Rib, 

E = ARig". (4) 

The exponent is usually found to lie in the range 1 < 
n < 1.75 and has been the subject of many theoretical 

studies of the entrainment process. Dimensional analy- 
sis suggests that Ris should be the controlling param- 
eter, so long as microscopic difl^usion plays a minor role 
(Phillips 1966). Basic energetic arguments in which the 
rate of change of potential energy due to mixing is as- 
sumed to be proportional to the turbulent kinetic energy 
available at the interface leads to an exponent of n = 1 
(e.g. Linden 1975). This same power law exponent has 



also been derived for models of the growth of the plan- 
etary boundary layer due to turbulent entrainment by 
penetrative convection (StuU 1973; Tennekes 1974; StuU 
1976a; Sorbjan 1996). 

The normalization of the entrainment coefiicient A has 
been found to vary significantly between the various lab- 
oratory and field studies conducted, with recent values 
found in the range 0.1 < A < 0.5 (e.g. Stevens & Brether- 
ton 1999). The discrepancy among the normalization 
constants has been called the 'A-dilema' (Bretherton et 
al. 1999). A review (up to 1991) of experimental mea- 
sures of the parameters in the entrainment law of equa- 
tion 4 arc tabulated in Fernando (1991) and a recent 
review of entrainment models used in the atmospheric 
sciences is discussed by Stevens (2002). 

The experimental and theoretical models discussed 
above arc generally motivated by geophysical problems, 
but are directly relevant to the conditions found in stellar 
interiors. The bulk Richardson numbers which charac- 
terize stellar convective boundaries fall within the same 
parameter range (10 < Ris < 500), and the background 
stratifications posses a similar buoyancy structure, so 
that it is interesting to learn from the geophysical models 
and compare to the stellar case. 

3. THE NUMERICAL TOOLS 

3.1. ID Stellar Evolution 

The hydrodynamic simulations which we study in this 
paper are of two distinct phases in the evolution of a 
23M0 supernova progenitor: main sequence core con- 
vection, and convective oxygen shell burning. The ini- 
tial conditions for the multi-dimensional simulations are 
taken from one-dimensional stellar models evolved with 
the TYCHO steUar evolution code. TYCHO (Young 
& Arnett 2005) is an an open source code^. A choice 
of standard ID stellar evolution procedures are used. 
The mixing length theory as described in Kippenhahn 
& Weigcrt (1990) is used with instantaneous mixing of 
composition in the convectively unstable regions. The 
limits of the convection zones are determined using the 
Ledoux criterion, which incorporates the stabilizing ef- 
fects of composition gradients. Semiconvective mixing 
has been turned off. Nuclear evolution is followed with 
a 177 clement network using the rates of (Rauscher & 
Thielemann 2000). Opacities are from (Iglesias & Rogers 
1996) and (Alexander & Ferguson 1994) for high and low 
temperature regimes, respectively. The solar abundance 
of GrevesHc & Sauval (1998) arc used. Although more 
recent abundance determinations have been made (As- 
plund et al. 2005) the impact on the stellar structure of 
the models presented here is small, and minor variations 
in the abundances have a negligible influence on the de- 
velopment of the hydrodynamic flow. 

3.2. Multi-Dimensional Reactive Hydrodynamics with 
PRO MP I 

The core of our multi-dimensional hydrodynamics code 
is the solver written by Fryxell, Miiller, & Arnett (1989) 
which is based on the direct Eulerian implementation 
of PPM (Colella & Woodward 1984) with generafiza- 
tion to non-ideal gas equation of state (Colella & Glaz 

^ http: //chandra. as . cirizona.edu/ dave / tycho-intro. html 
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1985). This code solves the Eiiler equations, to which we 
add nuclear reactions and radiative diffusion through an 
operator-split formulation. The complete set of combus- 
tive Euler equations, including diffusive radiative trans- 
fer, can be written in state- vector form, 



dQ 

'dt 

with the state vector 



Q 



+ v-* = s, 



p 

pu 
pE 
pXi 



(5) 



(6) 



the flux vector 



* = 



pn 

pVLU+p 

{pE + p)vi + 
pXiu 



(7) 



and the source vector 



S = 





/3U • g + penet 

Ri 



(8) 



where E ^ Ej + Ek is the total energy per gram con- 
sisting of internal and kinetic energy components, and p, 
p, u, g, and T are the density, pressure, velocity, grav- 
itational force field and temperature. The net energy 
source term due to nuclear reactions and neutrino cool- 
ing is Cnet = £burn + fcool, and the time rate of change 
of composition Xi due to nuclear reactions is denoted 
Ri. The radiative flux is = —krVT, with radiative 
"conductivity" = AacT^ j (Skup) and Rosseland mean 
opacity, k/j. Self gravity is implemented assuming the 
interior mass at each radius is distributed with spheri- 
cal symmetry. The mass interior to the inner boundary 
of the hydrodynamics grid is adopted from the TYCHO 
stellar model. 

The stellar models, which are calculated on a finely 
meshed Lagrangian grid, are linearly interpolated onto 
the Eulerian hydrodynamics grid taking into account the 
sub-grid representation of mass used in the PPM scheme. 
Mapping the models leads to small discrepancies in hy- 
drostatic equilibrium. An equilibration to hydrostatic 
balance occurs through the excitation and then damping 
of low amplitude, standing, predominantly radial pres- 
sure waves within the computational domain. These low 
amplitude waves, which are well described by the lin- 
earized wave equation, have a negligible affect on the 
convective flow. 

To save computational resources, we simulate care- 
fully chosen subregions of the star. Thus, these calcu- 
lations are local models of convection in the box in a 
star tradition. The advantage of local convection mod- 
els is that higher effective resolution can be used than 
is currently possible in global circulation models. This 
approach, however, precludes investigation of the lowest 
order modes of flow, and we do not yet include rotation 
or magnetic fields which are best studied using global 
domains. The boundary conditions used are periodic in 



angular directions, and stress- free reflecting in the radial 
direction. 

Our simulation code, dubbed PROMPI, has been 
adapted to parallel computing platforms using domain 
decomposition and the sharing of a three zone layer of 
boundary values and uses the MPI message passing li- 
brary to manage interprocess communication. 

4. OXYGEN SHELL BURNING 

We hiwc evolved a 23M0 stellar model with the TY- 
CHO code to a point where oxygen is burning in a shell 
which overlies a silicon-sulfur rich core. Approximately 
60% of the oxygen fuel available for fusion has been de- 
pleted at the time we begin the hydrodynamic simula- 
tion, when the star is ~ 2 x 10'' yrs from the zero age 
main sequence. Carbon, helium, and hydrogen burn- 
ing shells are also present contemporaneously at larger 
radii in the classic "onion skin" structure (Hoyle 1946). 
In one of the models presented here (ob.2d.e), which is 
also discussed in Meakin & Arnett (2006a), we adopt an 
outer radius that encompasses both the oxygen and car- 
bon burning shells. In this paper, however, we focus our 
analysis on the oxygen shell burning convection zone and 
the stable layers which bound it. 

The oxygen shell burning model affords us the oppor- 
tunity to study a thermally relaxed model because the 
thermal balance is determined by the very large neutrino 
cooling rates rather than the much lower radiative diffu- 
sion timescale (Arnett 1996, ch.ll). Neutrinos dominate 
the energy balance in the stable layers so that the stellar 
structure and the nature of convection are determined 
by the interplay between nuclear burning and neutrino 
emission (Aufderheide 1993; Arnett 1972). The eflPects 
of radiative diffusion are both unresolved and energeti- 
cally unimportant during these evolutionary phases, and 
has not been included in the oxygen shell calculations for 
computational efficiency. 

The radial profile of the simulated region is presented 
in Figure 2. Thc! temperature and density profiles be- 
tray the complex structure of the model, including the 
narrow burning shell that resides at the very base of the 
convection zone which is coincident with the tempera- 
ture peak. The initial extent of the convection zone 
can be identified by the plateau in oxygen mass frac- 
tion at 0.43 < rg < 0.72 (where rg = r/lO^cm). Char- 
acteristic of shell burning regions, the entropy gradient 
is quite steep at the boundaries of the convection zone 
and gives rise to peaks in the buoyancy frequency at 
those locations. The initial location of the upper con- 
vective boundary is coincident with a small stable layer 
at r ~ 0.72 x lO^cm, which is overwhelmed by the con- 
vective flow that develops in the simulation (see §4.1). 
A new boundary forms where the buoyancy frequency 
again becomes stabilizing at r > 0.8 x lO^cm. This mix- 
ing is shown in the change in ^^O abundance (Figure 2, 
top right) after 400s. 

In Table 1 we list the 25 nuclei used in our network. 
This network reproduces to within 1% the energy gener- 
ation of the full 177 element network used to evolve the 
one-dimensional TYCHO model for the simulated condi- 
tions, including oxygen and carbon burning shells. Dur- 
ing carbon burning the dominant reactions are ^^C(^^C, 
a)20Ne and ^^C(^^C, p)23Na, leaving an ash of ^^Ne, 
^^Na, protons and alpha-particles. ^"^Ne is photodis- 
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integrated through the ^''Ne(7, a)^^0 reaction. The 
dominant reactions during oxygen burning are ^^0(^®0, 
Q)28Si, 160(160, p)3ip, and 160(160, n)3iS, leaving an 
ash of predominantly ^^Si and ^^S. Neglecting the non- 
alpha chain species ^^Na, ^ip and ^iS can affect the net 
energy generation rate during carbon and oxygen burn- 
ing by a factor of a few under the conditions studied 
here. The reaction rates, including i^C(q!, 7)160, are 
from Rauscher & Thielemann (2000). 

Nuclear evolution is time advanced using the same re- 
action network subroutines as the TYCHO code and uses 
implicit differencing (Arnett 1996). We include cool- 
ing by neutrino- antineutrino pair emission, denoted ecooh 
which result from photo, pair, plasma, bremstrahlung, 
and recombination processes (Beaudet et al. 1967; Itoh 
et al. 1996). 

The Hclmholtz equation of state code of Timmes & 
Swesty (2000) is used to represent the ion and electron 
pressure with an arbitrary degree of electron degener- 
acy. With our 25 nuclei network, the initial conditions 
are thermodynamically consistent with the initial TY- 
CHO model to better than a few percent at all radii 
after mapping to the hydrodynamics grid. 

We calculate oxygen shell burning models in two and 
three dimensions. Our baseline model, labeled ob.2d.c, is 
a 90° wedge embedded in the equatorial plane with radii 
encompassing the oxygen burning convcctive shell and 
two stable bounding layers. The effects of dimensionality 
on the oxygen burning convective shell are explored with 
a three-dimensional model ob.Sd.B which has an angular 
extents of (27°)^. The influence of the upper boundary 
was studied with model ob.2d.e, which includes the over- 
lying carbon burning convective shell as well (additional 
details concerning this model are presented in Meakin 
& Arnett (2006a)). A preliminary resolution study is 
undertaken with model ob.2d.C which uses the same do- 
main limits but twice the linear resolution of the baseline 
model. Properties of the oxygen shell burning models 
presented in this paper are summarized in Table 2. 

4.1. The Correct Mixing Boundary 

Convection is initiated through random low-amplitude 
(0.1%) perturbations in density and temperature applied 
to a region in the center of the convectively unstable 
layer on a zone by zone basis. (Two additional simu- 
lation models with the same characteristics as ob.2d.c 
were calculated which used perturbations with a larger 
amplitudes and (1%), and a low order mode distribution. 
The development of the convective flow was found to be 
insensitive to these differences.) The role played by the 
perturbations is to break the angular symmetry of the 
initial model, and seed rising and sinking plumes whose 
growth is driven by nuclear burning, neutrino cooling, 
and the slightly superadiabatic background gradient im- 
printed in the initial TYCHO model. As the plumes rise 
they penetrate the original convective boundary which 
was determined in the TYCHO code using the Ledoux 
criterion. The initial evolution of the flow is presented 
in a time series of snapshots in Figure 3; the light yellow 
contour shows the initial outer convective boundary. 

The location of the initial outer boundary can be seen 
as a small bump in the initial profile of the buoyancy fre- 
quency presented in Figure 2 at radius r ~ .72 x lO^cm. 
The reason the boundary is stable in the ID model but 



did not survive in the multi-D simulation is because of 
the local nature of the Ledoux criterion used. This can be 
appreciated by the fact that although the buoyancy fre- 
quency at this location is positive, and hence locally sta- 
ble to convective turnover, the buoyancy jump across this 
region is very small Ab ~ 3 x 10^ cm/s^ compared to the 
turbulent kinetic energy in the adjacent flow, by which it 
is easily overwhelmed. This type of inconsistency can be 
relatively easily removed from ID simulations by using 
a parameter akin to the bulk Richardson number (eq. 
[1]) to characterize convective boundaries in place of the 
Ledoux or Schwarzchild criteria. For the original outer 
boundary Ris 1j ^ condition under which a boundary 
is expected to mix on an advection timescale, akin to the 
expansion of turbulence into a homogenous medium. 

The relationship between Ris and the traditional 
Schwarzschild and Ledoux criteria can be appreciated 
by writing the buoyancy frequency in terms of the well 
known "nablas" used in stellar evolution, 

iV^ = ;^(v„.-V. + |v,) (9) 

where V = (dlnT/dlnp), is the gradient of the stel- 
lar backgromid, Vad is the gradient due to an adiabatic 
displacement, = (rf In/i/dlnp) is the mean molecular 
weight gradient, and the thermodynamic derivatives are 
6 = —{dlnp/dlnT) and ip = (dlnp/dln/x). Therefore, 
the Ledoux criteria is simply, 

> 0. (10) 

The Schwarzschild criteria is the same, but with the sta- 
bilizing effect of the mean molecular weight gradient 
neglected. For comparison, the bulk Richardson number 
can be written Rib ^ N'^hL/a'^, where h is some mea- 
sure of the boundary width. A convective boundary will 
start to become stabilizing when, 

> a^/{hL). (11) 

This latter criteria is based on a finite threshold for sta- 
bility which takes into account the strength of the con- 
vective turbulence. Additionally, the bulk Richardson 
number is more than a simple stability criteria; it is also 
an indicator of the rate at which boundary erosion will 
proceed. We conclude that the correct criterion for de- 
termining the extent of a convective zone is neither the 
Ledoux nor the Schwarzschild criterion, which are both 
static, linear, and local criteria, but a dynamic boundary 
condition, based on the bulk Richardson number, which 
we will discuss in more detail in §7. 

4.2. Time Evolution 

The rich dynamics taking place at the convective 
boundary are apparent in the time evolution of the 3D 
flow presented in Figure 4, which provides a global view 
of the evolution. The upper panel shows the evolution in 
time and radius of the oxygen abundance gradient, rep- 
resented by a colormap in which light is large and dark 
is small. At the beginning of the simulation (far left) 
the colors are smooth as the turbulence has not yet de- 
veloped. The light line near the bottom of the panel is 
the lower boundary of the convective shell, where oxy- 
gen is separated from the silicon-sulfur core below. The 
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short horizontal band at r ^ 0.72 x 10 cm is the ini- 
tial weakly stable convcctive boundary discussed above; 
it is overwhelmed in the first 100 seconds by convection. 
After ~300 seconds the abundance distribution has ap- 
proached a quasi-steady state, with slow growth of the 
convective region. The bottom of the convection zone 
moves downward, but at a much slower rate than the 
upper boundary moves outward. The mottled apearance 
in the convection zone is due to the ingestion of new oxy- 
gen entrained from above, followed by turbulent mixing. 
At the top boundary of the convection zone an oscilla- 
tory behavior can be seen, and in the overlying stable 
region wave motions arc apparent. 

The lower panel in Figure 4 shows the radial profile of 
the kinetic energy, which illustrates a major feature of 
the convection: intermittency. While these simulations 
are well described by a statistical steady state over a few 
convective turnover times, at any instant the fluctuations 
are significant. The flow is episodic, with bursts of activ- 
ity followed by lulls. The bursts in kinetic energy in the 
convection zone are seen to induce wave trains in both 
the upper and lower stable layers. Characteristic of g- 
modes, the phase velocity (orientation of the wave crests) 
is orthogonal to the group velocity (direction of energy 
transport) in these wave trains, which can be seen by 
comparing the composition and kinetic energy profiles. 

4.3. Quasi-Steady Oxygen Shell Burning Convection 

Following the transient readjustment of the outer con- 
vective boundary, the oxygen burning convective shell at- 
tains a quasi-steady character. In Figure 5 we present the 
time evolution of the integrated internal, gravitational, 
and kinetic energy. The energy is calculated by form- 
ing horizontal averages of the flow properties and then 
assuming a full spherical geometry. The gravitational 
energy contribution from material on the computational 
grid is calculated according to, 

Ea^J^^^dr, (12) 

where the mass increment is dM = 47rr^(/9), and the 
integral is taken over the radial limits of the grid. 

The total kinetic energy levels off in all of the mod- 
els by t ~ 300s. The 2D models are characterized by a 
much larger overall kinetic energy. The total kinetic en- 
ergy settles down to a slow increase as the oxygen shell 
evolves; this is true for both 2D and 3D. 

The radial profile of the r.m.s. velocity fluctuations are 
presented in Figure 6 for the 2D and 3D models. The ve- 
locity fluctuation amplitudes in all of the 2D models are 
higher than the 3D model by a factor of ^2. The 2D 
models also assume a significantly different radial profile 
than the 3D model, with a flow structure that is domi- 
nated by large convective vortices which span the depth 
of the convection zone. The signature of these large ed- 
dies is apparent in the horizontal velocity components, as 
well as the fairly symmetric shape of the radial velocity 
profile within the convection zone. The velocity c;ompo- 
nents in the 3D model reveal an up and down flowing 
circulation with horizontal deflection taking place in a 
fairly narrow layer at the convective boundaries. 

Although significant differences exist between 2D and 
3D models, the 2D models are found to be in good agree- 



ment with each other to the extent that the statistics 
have converged, which are calculated over the time pe- 
riod t G [300,450]s. The time period for calculating 
statistics was limited by the model ob.2d.C, which was 
only run as far as f ~ 450s. The agreement among 
the 2D models shows that the outer boundary condition 
(tested by model ob.2d.e) and the grid resolution (tested 
by model ob.2d.C) are not playing a decisive role in de- 
termining the overall structure of the flow, at least in 
these preliminary tests. The agreement in overall veloc- 
ity amplitude in the upper stable layer in model ob.2d.e 
indicates that the stable layer velocity amplitudes arc 
not strongly affected by the details of the modes that are 
excited in that region. This gives credence to the analy- 
sis in (Meakin & Arnett 2006b) which assumes that the 
stable layer velocity amplitudes are determined by the 
dynamical balance between the convective ram pressure 
and the wave induced fluctuations. 

The convective turnover times tc = 2AR/vconv for the 
2D models are all of order tc ^ 40s, and they span be- 
tween 10 and 55 convective turnovers. The turnover time 
for the 3D model is tc ^ 103s, and the model spans ap- 
proximately 8 convective turnovers. 

4.4. Stable Layer Dynamics During Shell Burning 

In both the 2D and 3D models, the stably stratifled 
layers are characterized by velocity fluctuations through- 
out their extents (Figure 6). These fluctuations are the 
signature of g-modes which are excited by the convective 
motions. In the 2D model, the amplitudes of the stable 
layer velocity fluctuations are higher. In the lower stable 
layer, the 2D models also have a much smaller ratio of 
horizontal to radial velocity amplitude. The velocity am- 
plitude ratio is roughly proportional to the ratio of the 
mode frequency and buoyancy frequency, Vr/v± w w/N 
(Press 1981), so that the waves excited in the 2D model 
are of lower frequency. The velocity ratios in the upper 
stable layer are comparable between the 2D and 3D mod- 
els, though the 2D amplitudes are higher by a factor of 
~ 2. 

During late burning stages, multiple concentric convec- 
tive shells form which are separated by stably stratified 
layers. These intervening stable layers act as resonating 
cavities for g-modes that are excited by the turbulent 
convection. In Meakin & Arnett (2006a) it was shown 
that the stable layer motions in model ob.2d.e can be de- 
composed into individual g-modes that are well described 
by the linearized non-radial wave equation (Unno et al. 
1989). Meakin & Arnett (2006b) showed that a good es- 
timate for the amplitudes of the wave motions (and the 
associated thermodynamic fluctuations) in both the 2D 
and 3D models can be made by assuming that the pres- 
sure fluctuations associated with the g-modes balance the 
ram pressure of the turbulent convection. In the latter 
paper, a single mode (frequency and horizontal scale) was 
assumed, based on integral properties of the turbulence 
(convective turnover time, and mixing length scale). In 
this section we present the spectrum of motions present 
in the stable layers and turbulent regions for the more 
realistic 3D model. 

For a given background structure, a spectrum of eigen- 
modes exist which are solutions to the non-radial wave 
equation and boundary conditions. Individual modes can 
be uniquely identified by a horizontal wavenumber index 
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/ and oscillation frequency lo. In Figure 7, I — lu dia- 
grams arc presented for the convection zone, and the two 
bounding stable layers. The individual / — to components 
have been isolated through Fourier transforms of a time 
sequence of the simulation data. 

Several modal components or "branches" can be iden- 
tified in the stable layer diagrams (left and right panels). 
These include: (1) p-modes, seen as a series of points at 
the lowest I values that extend to high frequencies; (2) g- 
modes, which appear as ridges that are bound above by 
the buoyancy frequency; and (3) f-modes, which appears 
as a ridge separating the g- and p-modes. The /-modes 
are interfacial waves, and are most prominently seen in 
the lower-boundary diagram at a radius r = 0.4 X 10® 
cm. The f-mode signature is due to interfacial waves run- 
ning along the convective boundary at r ~ 0.43 x lO^cm, 
where there is a spike in buoyancy frequency. 

In the convection zone, the spectrum is dominated by 
power at low temporal and spatial frequencies. This 
strong non-modal convection signature is also present, 
though at lower amplitude, in the stable layers. This 
"turbulence" spectnun can be seen extending from the 
lower left corner of the diagrams. This same feature was 
also present in the simulations of He-shell burning by 
Herwig et al. (2006). 

5. CORE CONVECTION 

Are the hydrodynamic features of oxygen shell burn- 
ing of more general applicability? To investigate this, 
we examined core convection during hydrogen burning. 
Because of the long thermal time scale for radiative dif- 
fusion in such stars, we focus on the hydrodynamic be- 
havior of a model in which the inner boundary provides 
a driving luminosity about ten times larger than natural. 
This allows us to simulate the convection with our com- 
pressible hydrodynamics code; an anelastic method (if 
multi-fluid) would allow this to be done in the star's nat- 
ural time scale. While our calculation is not thermally 
relaxed on a Helmholtz-Kelvin time scale, it does relax 
dynamically, and provides some clue as to the convective 
behavior. 

We have previously evolved a 23 M© star onto the 
main sequence with TYCHO, to an age of 2.4x10^ yrs, 
at which point hydrogen is burning in a convective core. 
The model is then mapped onto the PROMPI hydro- 
dynamics grid for simulation. This model represents an 
early point in main sequence evolution, in which the core 
hydrogen content has been depleted by only 1.7% {Xcore 
= 0.689, Xinit = 0.701, AX = 0.012). The inner radius 
of the simulation was chosen such that the convectively 
unstable region covers ~1 pressure scale height (convec- 
tive cores are usually only of order a pressure scale height 
because of the divergence of the scale height towards the 
stellar center). The entire domain covers ~5 pressure 
scale heights and '~3.3 density scale heights. The density 
contrast across the computational domain is ~30 with a 
contrast of ^2 across the convectively unstable region. 

The radial profile of the simulated region is presented 
in Figure 8 including the run of temperature, density, 
composition, buoyancy frequency, and relative buoyancy. 
The entropy jump at the edge of the convective core, due 
to the fuel-ash separation, gives rise to a buoyancy jump 
(spike in buoyancy frequency). 

The equation of state for the main sequence model is 



well described by an ideal gas with radiation pressure 
component. The ratio of gas to total pressure lies in 
the range 0.85 < /3 < 0.95, with an increasing contribu- 
tion from radiation pressure as the stellar center is ap- 
proach. A single composition representing hydrogen has 
been evolved to keep track of nuclear transmutation and 
the mean molecular weight of the plasma. A metalliicity 
of Z — 0.01879 has been used to represent the additional 
175 species in the initial TYCHO model and helium is 
calculated according ioY = 1— [X -\- Z), where X is the 
self consistently evolved hydrogen mass fraction. 

The luminosity due to nuclear burning in the compu- 
tational domain is a small fraction of the total stellar 
luminosity (2.4%) which is dominated by burning in the 
inner regions of the core and Ltot = 7.8 x IQ'^Lq. Core 
burning is incorporated into the simulation as an input 
luminosity at the inner boundary of the computational 
domain. 

The Kelvin-Helmholtz timescale for this model is 
tKH ~ 10^ years, which is many orders of magnitude 
longer than the dynamical timescales that are feasible to 
simiilate. Additionally, the small luminosity of the star 
produces a convective velocity scale that is very subsonic 
(M ~ 10"'^). Since we are not interested in the thermal 
relaxation of the model, we have boosted the input lu- 
minosity by a factor of 10 to increase the velocity scale 
of the flow. This was necessary because our fully com- 
pressible code is limited by the sound crossing time. (An 
anelastic or low-mach number method would be ideal for 
simulating this core convection flow at the natural veloc- 
ity scale.) 

Radiation transport is treated in the diffusion limit. 
Opacities approximated by Thompson scattering, which 
agrees well with the OPAL opacities (Iglesias & Rogers 
1996) used in the ID TYCHO model for the region sim- 
ulated. The efli'ects of radiative diffusion, however, are 
found to be unresolved in the current simulation (with 
the diffusion time across a single zone Trad — /K-ad ^ 
tconv, with grid zone size A) and therefore energy trans- 
port in the convection zone occurs primarily on the sub- 
grid level due to numerical diffusion. This is a high Peclet 
number simulation. 

A 2D and a 3D model have been calculated. The sim- 
ulated wedges have angular extents of 30° in both the 
polar and azimuthal directions and are centered on the 
equator to avoid zone convergence problems near the 
poles. This minimal angular domain size was chosen 
by calculating models of increasing angular size in 2D 
domains until the flow pattern converged. The angular 
domain size used in the present simulations encompasses 
a large convective roll in 2D. Smaller 2D domains were 
found to distort the convective roll while domains larger 
by integer multiples contained proportionally more rolls 
of the same flow amplitude and morphology. The bound- 
ary conditions in the radial direction are reflecting and 
stress free, and periodic conditions are used in both an- 
gular directions. The grid zoning, domain limits, and 
simulation run times are summarized in Table 3 for the 
2D and 3D models. 

5.1. Quasi- Steady Core Convection 

Convection is initiated through random low-amplitude 
(0.1%) perturbations in density and temperature applied 
as in the oxygen shell simulation. In order to save com- 
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puting time, the 3D model was initiated on a domain one 
quarter as large in azimutlial angle which was then tiled 
four times in angle once convective plumes began to form. 
The initial development of the flow in the 3D model is 
presented in Figure 9 as a time sequence of velocity iso- 
surfaces. The turbulent structure of the convective flow, 
as well as the excitation of internal waves which radiate 
into the overlying stably stratifled layer, are clearly illus- 
trated. A comparison of the flow morphology between 
the 2D and 3D models is presented in Figure 10. The 2D 
convective flow is much more organized and laminar, and 
is dominated by a single large convective cell while the 
3D convection is composed of many smaller scale plume- 
like structures and is more obviously turbulent. In both 
models the stably stratified regions are rife with internal 
waves excited by the convection. 

The 3D convective flow attains a quasi-steady charac- 
ter after approximately 6x10^ s, or approximately two 
convective turnovers. The evolution of the internal, grav- 
itational, and total kinetic energy components on the 
computational grid for the 2D and 3D models, are pre- 
sented in Figure 11 and are calculated in the same way 
as for the oxygen burning model. 

In both the 2D and 3D models, the total kinetic en- 
ergy fluctuates in times with excursions from the mean as 
larger SEk/^ ^ 0.4 in 3D and 5Ek/'E^ ^ 0.6 in 2D. 
The kinetic energy in the 2D model grows on a slightly 
longer timescale, and achieves a steady character after 
^ 10^ s, at which time the kinetic energy growth rate 
tapers off. The total energy is conserved to better than 
~0.2% for both the 2D and 3D flows by the end of the 
calculation. 

The r.m.s. velocity fluctuations are presented in Fig- 
ure 12 for the 2D and 3D models. The resultant flows 
in both the 2D and 3D models are similar to that found 
for the oxygen shell burning model. The velocity ampli- 
tudes are higher in 2D by a factor of ^5 (sec axis scale 
in Figure 12), and the flows are dominated by large ed- 
dies spanning the depth of the convective region. The 
horizontal deflection of matter is also found to occur in 
a much narrower region in the 3D model. The hard-wall 
lower boundary of the 3D model is characterized by an 
even narrow horizontal flow, probably due to the absence 
of a stable layer which is host to g-modes. 

The time averaged convective flow velocity for the 3D 
model is Vc ~ 2.8 x 10^ cm/s. The turnover time is 
tc = 2AR/vc ~ 3.2 X 10''' s, and the simulation spans ap- 
proximately 5 convective turnovers. The peak velocity 
fluctuation is Vpeak ~ 2 X 10^ cm/s, corresponding to a 
peak Mach number of M ~ 0.03, and the maximum den- 
sity fluctuations within the convective flow are ~ 0.02%, 
which is of order as expected for low Mach number 
thermal convection (Gough 1969). The time averaged 
convective flow velocity in the 2D model is ~ 1.3 x 10^ 
cm/s, and the convective turnover time for this model is 

~ 7 X 10^ s. The simulation spans 1.5 x 10^ s which is 
^21 convective turnovers. The peak velocity fluctuation 
in the 2D model is comparable to that in the 3D simu- 
lation, with Vpeak ^ 2 X lO*' cm/s and the peak density 
fluctuations is a little more than twice that found in the 
3D model, -^0.05%. The turnover times and convective 
velocity scales are summarized in Table 3. 



5.2. The Stable Layer Dynamics Overlying the 
Convective Core 

As in the oxygen shell burning model, the stably strat- 
ified layers in the core convection models arc charac- 
terized by velocity fluctuations throughout. Similar to 
shell burning, the 2D stable layer velocity amplitudes are 
larger and have a smaller radial to horizontal component 
ratio Vr/v± « u/N compared to the 3D flow. 

The stable layer motions in the core convection simula- 
tion arc predominantly resonant modes, which compare 
well to the analytic eigenmodes of the linearized wave 
equation, and are analogous to those discussed for the 
oxygen shell burning model. The region outside the con- 
vective core will act as a resonant cavity, with the outer 
boundary at the location where the buoyancy frequency 
and Lamb frequency cross. 

The amplitudes of the internal waves excited will 

be determined by the ram pressure of the turbulence at 
the convective boundary. In Figure 13 the ram pressure 
and horizontal r.m.s. pressure fluctuations are presented 
for the 3D model, and can be seen to balance at the inter- 
face between the convective core and the stably stratified 
layer. Using this condition of pressure balance, Meakin 
& Arnett (2006b) estimate the amplitudes of the excited 
internal wave velocities and the induced thermodynamic 
fluctuations and find this to be in good agreement with 
the oxygen shell burning simulations. The relationship 
between the density fluctuations, the convective velocity 
scale, and the stellar structure (i.e., N and (?) was found 
to be, 



That this proportionality holds in the core convection 
model as well, where fluctuation amplitudes are lower 

than those in the oxygen shell burning model by an order 
of magnitude, is illustrated in Figure 14, which presents 
the buoyancy frequency and density fluctuation profiles 
for the boundary region. The measured density fluctu- 
ation and the value calculated according to equation 13 
compare remarkably well, with p' /{p) ^ 0.12%. 

6. SIMULATIONS AND MIXING LENGTH THEORY 

In this section we compare our 3D oxygen shell burning 
simulation results to the mixing length theory of convec- 
tion. We choose to compare this particular simulation 
since it represents the most physically complete model 
in our suite of calculations, both in terms of dimension- 
ality and thermal evolution. Unless otherwise specified, 
the time period over which averaging is performed on 
the simulation data is t £ [400, 800] s, which is approxi- 
mately 4 convective turnovers. We find that this period 
is sufficiently long compared to the time evolution of the 
flow that average values are not affected appreciably by 
increases in the averaging time period. 

6.1. Mixing Length Theory Picture 

The basic picture underlying the mixing length the- 
ory, which is the standard treatment of convection used 
in one-dimensional stellar evolution modeling (see Cox & 
Guifi 1968; Clayton 1983; Kippenhahn & Weigert 1990; 
Hansen & Kawaler 1994), is one in which large eddies 
are accelerated by an unstable temperature gradient and 
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advect across a certain distance until they suddenly lose 
their identity by turbulently mixing with the background. 
Energy is transported through this process because the 
envisioned turbulent blobs which are moving radially out- 
ward have a higher entropy at their formation location 
than the location in which they dissolve. The vertical ex- 
tent over which large eddies retain their identity as they 
advect through a convection zone is a fundamental pa- 
rameter in the mixing length theory. This mixing length 
A is generally taken to be a multiple of the local pressure 
scale height K = at^Hp. 

Within this physical picture, the mixing length the- 
ory develops a relationship between the convective flux, 
the temperature gradient, the velocity scale, and the ge- 
ometrical factors which describe the large scale eddies. 
The starting point in mixing length theory is the radial 
enthalpy flux, which is written in terms of fluctuations 
in the flow properties, and is taken to be (assuming a 
horizontally isobaric flow), 

F, ^ VcpcpT'. (14) 

The temperature fluctuations in mixing length theory 
are related to the temperature gradient and the distance 
traveled by. 



important additional parameter since the eddies are en- 
visioned to leak a fraction of their thermal energy over a 
mixing length distance. When local cooling dominates, 
either through radiative losses (in optically thin regions) 
or neutrino losses (such as in the present model), the 
geometry of the eddies is not important since energy es- 
capes everywhere from the large eddies, not just at eddy 
" surfaces" . 

During the oxygen shell burning simulations being con- 
sidered here, non-adiabatic losses are small over a con- 
vective turnover time and the convection is expected to 

be "efficient". A quantitative measure of convective ef- 
ficiency is the Peclet number, which is the ratio of the 
energy loss timescale to the convective turnover timescale 
for the large eddies. In the current model, energy losses 
are dominated by neutrino cooling e,^. Therefore, we 
calculate an effective Peclet number using the following 
convective and neutrino-cooling timescales: 



" T'de^/dT ~ KdlaTJ ^ ' 



T'/T = 



dhiT, dlnTnxA 



dr 



dr 



(AV) 



1 A 

Hp 2 



(15) 



where the subscript "e" indicates properties of the large 
convective eddies and the dimensionless temperature gra- 
dient V is used (see §4.1), and the difference between the 
gradient in the eddy as it moves and the averaged stellar 
background is written, 

AV = (V - Ve). 

The factor of 1/2 in equation 15 represents the idea that 

on average the large convective eddies have traversed 
about half a mixing length before reaching the current 
position. The velocity obtained by the convective eddy is 
computed by calculating the work done by the buoyancy 
force over a mixing length. 



vl = ff/3(AV) 



A^ 

8Hr, 



(16) 



Here again, the eddy is assumed to have been accelerated 
over half of a mixing length and an additional factor of 
1/2 is incorporated on the right hand side to account 
for energy lost driving other flows, such as small scale 
turbulence and horizontal motions (e.g., note that the 
r.m.s. horizontal velocity is of the same order as the 
r.m.s. radial velocity in the simulation). The average 
convective flux can then be written. 



= pCpT^/^ 



A2 



(AV) 



3/2 



(17) 



The temperature gradient for the convecting material 
is found by assuming that eddies follow isentropic tra- 
jectories Ve = Vad- Deviations from isentropic motion 
have been considered in the mixing length theory. In the 
case of strong radiative diffusion losses, the eddy geom- 
etry (in terms of the surface area to volume ratio) is an 



Pe = 



VcCpT /(?lne,y 



dlnT 



Ftpf'v 



where characteristic values from the simulation have 
been used in equation 20, and the temperature depen- 
dence of the neutrino loss rates is (91ney/91nT) < 9. 
Therefore, the Peclet number for the convection is Pe 
> 10^, and we should expect the convection zone to be 
very nearly isentropic. 

6.2. The Enthalpy Flux, Background Stratification, and 
Temperature and Velocity Fluctuations 

The convective enthalpy flux measured in the simula- 
tion is presented in Figure 15. The spike at the bottom 
of the convection region and the slight dip at the top 
reflect the braking of convective motion at these bound- 
aries. The enthalpy flux is calculated by performing time 
and horizontal averages on the flow. Mixing length the- 
ory, however, makes the assumption that the velocity 
and temperature fluctuations are perfectly correlated, so 
that horizontal averaging of fluctuations is comparable 
to products of the averages. This is not neccesarily the 
case. To test the degree to which the velocity and tem- 
perature fluctuations are correlated we can calculate the 
correlation coefficient, Ue, deflned in the following way: 



F, = {pCpT'v',) = aE{pCp){T'){v',). 



(21) 



The fluctuations, T' and v'^ in equation 21 are taken to 
the be the r.m.s. fluctuations in the simulation. The 

radial profile of ue is shown in Figure 15. We find 
{ue) = 0.7 ± 0.03 averaged over the radial interval 
r e [0.5, 0.75] X 10'^ cm within the convection zone. A 
value smaller than unity indicates that the horizontal dis- 
tribution of temperature and velocity fluctuations in the 
flow are not perfectly correlated. The degree of correla- 
tion, however, is found to be fairly uniform throughout 
the convection zone. 
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We next consider how well the velocity and temper- 
ature fluctuations are correlated with the local temper- 
ature gradient. In Figure 16 the temperature gradient 
of the horizontally averaged hydrodynamic model profile 
Vs as well as the adiabatic Vad and the composition- 
corrected (Ledoux) gradient V Led = ^ad + <f/P ^ V^, 
are presented. The super-adiabatic temperature profile 
of the stellar background AV = Vs — Vad is presented in 
the right panel of Figure 16. While the convection zone is 
found to have a super-adiabatic profile throughout, it is 
very small (AV < 10""^). This confirms the efficiency of 
the convection, in accord with our estimate for Pe. Sta- 
bility is maintained in the upper boundary layer by the 
composition gradient, where we have Vad *^ V^ <C ^ Led- 

In order to assess the validity of the mixing length 
theory temperature and velocity fluctuation amplitudes 
given by equations 15 and 16 we calculate the correlation 
coefficients ax and a„ which are defined by, 

T'/T = (AV)aT (22) 

and 

v, = ^^g(i{AW)H,. (23) 

An important question concerns how to interpret and 
measure the temperature and velocity fluctuations T' 
and Vc in the simulations for comparison to mixing length 
theory. In the mixing length theory, these fluctuations 
are identified with the properties of large eddies. There- 
fore, a direct comparison would entail isolating the large 
eddies from the rest of the flow, and measuring their 
properties. In lieu of this complicated procedure we iden- 
tify the fluctuations in the large eddies with two distinct 
quantities for comparison: (1) the r.m.s. fluctuations in 
the flow; (2) the difference between the horizontally av- 
eraged background value and the mean values in the up 
and down flowing material. 

The temperature fluctuations calculated using these 
two methods are presented in Figure 17. The temper- 
ature fluctuations in the convection zone follow a trend 
similar to the super-adiabatic gradient, i.e., decreasing 
with increasing radius. In the right panel of Figure 17 
the radial proflle of ax is shown using both deflnitions of 
the fluctuations. The nonzero temperature fluctuations 
outside the convective region are due to distortions in sta- 
ble layers due to convective buoyancy braking Meakin & 
Arnett (2006b) ; the use of separate up and down flows is 
cleaner, eliminating these. The slope in the temperature 
fluctuation proflles are slightly overcompensated for by 
the super-adiabatic gradient when forming the ratio ar- 
Within the scatter, however, ut is fairly well represented 
by a constant value. The mean value within the body of 
the convection zone (taken to be r e [0.5, 0.75] x 10^ cm) 
is larger for the r.m.s. fluctuations (aT(i'nis)) = 0.73 
compared to (aT(up)) = 0.45 and (aT(down)) = 0.40. 
The largest departures from the mean, within the con- 
vective region, occur at the base of the convection zone, 
r < 0.52 X 10^ cm, where the rniclear flame is driving 
the convective flow. The departures are also large at the 
top, in the region of buoyancy braking. 

The corresponding analysis for the velocity fluctua- 
tions is presented in Figure 18. The overall trends 
are similar for ar and a„. The mean values of ay 



within the body of the convection zone are (Q:-„(rms)) = 
1.22,(a„(up)) = 1.08, and (a„(down)) = 0.96. The 
largest departure from constancy is again found to be 
at the base of the convection zone (the flame region). 

The sharp decrease in the effective mixing length near 
the lower boundary is not entirely surprising. The dis- 
tance to the convective boimdary provides an upper limit 
to the mixing length, while further away from the bound- 
aries the mixing length is limited by the distance over 
which eddies can maintain their coherence. This effect 
is possibly more exaggerated at the lower boundary be- 
cause of the steep gradient in velocity which is needed to 
move the energy out of the burning zone. In contrast, the 
upper boundary is characterized by a more gentle decel- 
eration of material and a "softer" boundary (i.e., lower 
A''^). Ignoring this boundary effect and using the same 
mixing length parameter throughout the convection zone 
would result in a shallower temperature gradient near the 
boundary. The stiff temperature dependance of the nu- 
clear reaction rates may therefore be affected. 

The absolute calibration of ar and a„ are somewhat 
arbitrary, and are scaled by factors of order unity for 
a particular implementation of the mixing length the- 
ory based on the heuristic arguments discussed above. 
According to equations 15, 16, 22, and 23 the equiv- 
alencies are Q!a,t = 2 x cct and q;a,i; = x a„ 
where the values subscripted by A indicate the mix- 
ing length theory values defined by Kippcnhahn & 
Weigert (1990). The corresponding values measured in 
the simulation are (aA,v(rms)) = 1.73, (q;a_^(up)) = 
1.53, and (Q!A,t,(down)) = 1.35, for velocity fiuctua- 
tions; and (aA,T(rms)) = 1.46, (aA,T(up)) = 0.9, and 
(Q:A,T(down)) = 0.8 for temperature fiuctuations. 

The ratios Q!A,T/aA,ij are 0.84, 0.59, and 0.60 for the 
r.m.s., up-flow, and down-flow values, respectively. In 
relation to the present simulation, a higher degree of 
consistency (i.e., aA.u = Q^a.t) can be brought to this 
implementation of the mixing length theory by scaling 
the velocity fluctuation in equation 16 by the inverse of 
the ratio aA,T/oiA,v Physically, this translates into a 
higher efficiency (by a factor ~1.2 - 1.7) for the buoy- 
ancy work to accelerate the large eddies over the value 
of 1/2 adopted above, which is reasonable considering 
the heuristic argument used. Alternatively, agreement 
can be made by scaling the temperature fiuctuations in 
equation 15 by the same ratio, which amounts to decreas- 
ing the distance over which eddies remain coherent and 
adiabatic as they move across the convection zone. Both 
of the these effects are plausible, as well as a combination 
of the two so long as the ratio is maintained. Which is 
operating in the present simulation? Unfortunately, the 
degeneracy between these two parameters cannot be bro- 
ken because they combine linearly when calculating the 
enthalpy flux, which therefore docs not provide a further 
constraint. Finally, it is possible that the effective mix- 
ing lengths for temperature and velocity fluctuations are 
different, a notion that is supported by the correlation 
lengths which we discuss next. 

6.3. Correlation Length Scales 

In the top two panels of Figure 19 the vertical cor- 
relation length scales, calculated according to equation 
Bl, are presented for the velocity and temperature fluc- 
tuations. The vertical scale height is defined as the full 
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width at half maximum of the correlation function, and 
can be written in terms of the correlation length in the 
positive and negative directions, Ly = Ly — Ly. The 
relative values of Ly and Ly give an indication of asym- 
metries in the eddies (Figure 19, lower-left): L'^r/Ly = 1 
is a symmetric eddy; Ly/Ly > 1 is an eddy flattened 
on the bottom; and Ly/Ly < 1 is an eddy flattened 
at the top. Based on this simple diagnostic both the 
temperature and velocity correlations indicate that the 
eddies near the lower boundary are flattened on the bot- 
tom, and those at the upper boundary are flattened on 
the top. The "overshooting" distance (h ~ 10^ cm at 
the upper boundary and h < 10*' cm at the lower bound- 
ary), which is best described as an elastic response to the 
incoming turbulent elements, is very small compared to 
the correlation lengths measured here. Therefore, these 
eddies are effectively hitting a "hard wall" upon reaching 
the boundaries. 

The signature of this eddy " flattening" is also present 
in the radial profile of the full width length scale, Ly- In 
the case of velocity, which has larger correlation length 
scales, significant asymmetries arc present throughout 
the convection zone. The smaller length scales associ- 
ated with the temperature fluctuations permit a broad 
region throughout the convection zone where the eddies 
are roughly symmetric {Ly/Ly « 1) and appear to be 
uninfluenced by the boundaries. In this intermediate re- 
gion, away from the boundaries, the temperature fluctu- 
ation length scales are relatively constant in size, even 
decreasing with radius, in contrast to the pressure and 
density scale heights which are increasing with radius. 

In the standard mixing length theory, the sizes of con- 
vectivc eddies arc assumed to be comparable to the size 
of the mixing length. How do the correlation length 
scales compare to the mixing length parameters found 
above? The lower left panel of Figure 19 shows the 
ratios of Ly to the pressure and density scale heights. 
None of these curves are particularly constant within 
the convection zone, and boundary effects are partic- 
ularly strong throughout the convection zone in the 
case of the velocity correlations. Interestingly, the ve- 
locity correlation parameter av{vr,Hp) = Lv{vr)/Hp 
is larger than the temperature correlation parameter 
a^{T',H.p) = Lv{T')/Hp. This is in accord with the 
ratio a\:r / (y:\,v < 1 found in the mixing length analy- 
sis above. Concerning the absolute calibration, however, 
the correlation length scales are smaller than the mixing 
length values by as much as a factor of a few. In an anal- 
ogous comparison by Robinson et al. (2004) for subgiant 
atmosphere models, the vertical correlation lengths were 
also found to be smaller than the mixing length used 
to construct the initial model, and that the ratio varied 
significantly throughout the convection zone. 

The horizontal correlation lengths Lh are shown in 
Figure 20 together with the vertical scales for compari- 
son. For the velocity, the horizontal scale is much smaller 
than the vertical, indicative of eddies which are signifi- 
cantly elongated in the vertical direction. The temper- 
ature fluctuations appear to be much more symmetric, 
with only a small degree of elongation in the vertical di- 
rection which is slightly more pronounced near the top of 
the convection zone. In the stable layers, the horizontal 
scales are larger than the vertical, which is a characteris- 



tic of the horizontal " sloshing" motions associated with 
g-modes. 

6.4. The Kinetic Energy Flux, Flow Asymmetry, and 
Moving Beyond the Mixing Length Theory 

The kinetic energy flux associated with convection is 
ignored in the mixing length theory since it arises from 
the asymmetries in the flow and MLT assumes that the 
flow is symmetric. An order of magnitude estimate for 
the kinetic energy flux, however, can be made: 



Fk pvl/2 Vc 



Fc pCpT' Vc 



8 Tpcp 8 



^Vad ~ 0.03 (24) 



where mixing length relationships have been used to cal- 
culate Vc and T', a\ is assumed to be of order unity, and 
Vad ~ 0.25 has been adopted from the simulation. This 
result tells us that the kinetic energy flux will be a few 
percent of the convective enthalpy flux. This estimate is 
an upper limit since up-flows and down-flows will cancel 
to some degree (Bohm-Vitense 1992, §6.1). In the sim- 
ulation, the ratio of kinetic to enthalpy flux is found to 
be Fx^max/ Fc,max ^ 0.01, which is of order the simple 
MLT scaling, but down by a factor of a few as expected. 

We can directly relate the kinetic energy flux to the 
flow asymmetry in the following way. The upflow area 
covering fraction /„ = A^p/Atot is shown in Figure 21. 
We can then write an estimate for the kinetic energy flux 
as, 



FK,net = -^PoifuVl " fd^d) 



(25) 



which can be written in terms of just the flow velocities, 



1 



K.net 



Vu/Vd + 1 



v'd 



(26) 



where we have used the mass conservation equation, 
fuVu + fdVd = assuming p„ « pd which is a good 
approximation in these simulations. The kinetic energy 
flux in the simulation is shown in Figure 22. Shown by 
the thin line is Fk calculated according to equation 26, 
which is in good agreement. Here, we have used the hor- 
izontal and time averaged values for (v) and (v^). The 
mixing length theory, however, does not provide infor- 
mation about {v^)j but only (v). We overplot with the 
dashed line Fk calculated using {v)^ in place of (w^); it 
has to be scaled by a factor of 5 to fit the simulation 
data. 

The scaling factor needed to calculate the kinetic en- 
ergy flux is required to account for the skewness in the 
radial velocity field. More precisely, the correlation coef- 

3 

ficient X = is needed, which is related to the 

skewness 7 = (f^)/cr^. Both x and 7 are presented in 
Figure 22. Note, that the skewness is a good proxy for 
the down-flow covering fraction (/^ = 1 — /«; see Figure 
21), and therefore its sign is indicative of the direction of 
the kinetic energy flux. 

Convective regions which are spanned by several pres- 
sure scale heights are found to have kinetic energy to 
enthalpy flux ratios larger than a few percent as found in 
this study. For instance, the simulations of Cattaneo et 
al. (1991) and Chan & Sofia (1989) which each span -^5 
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pressure scale heights achieve \Fk / Fc\ ^ 35%, and the 
domain in Chan & Sofia (1996) spanning ~ 7 pressure 
scale heights achieves \Fk/Fc\ ~ 50%. A key result in 
the analysis of Cattaneo et al. (1991) is that the kinetic 
energy flux is dominated by coherent, downward-directed 
flows which are correlated over distances comparable to 
the simulation domain. Additionally, the enthalpy flux 
and kinetic energy fluxes associated with these down- 
flows essentially cancel with CpT' ^ v"^ , which was shown 
to follow if the downflows can be described as Bernoulli 
streamlines. 

The long range correlations just described, together 
with the boundary effects which dominate our shell burn- 
ing model, undermine the basic mixing length theory pic- 
ture of convection. The large coherence of the flows seen 
in these simulations, however, and present even in tur- 
bulent parameter regimes, suggest that modeling these 
coherent structures is a viable approach. Already, mod- 
els incorporating multiple streams or " plumes" as closure 
models (e.g. Rempel 2004; Lesaffre et al. 2005; Belkacem 
et al. 2006) are providing enticing alternatives to the mix- 
ing length theory. 

6.5. Related Studies 

Although the mixing length parameters calculated 
above deviate from constancy near the convective bound- 
aries, a mean value is a good approximation for most of 
the convection zone. It would be interesting if these pa- 
rameters a-i\ and ciy were universal, as assumed by 
mixing-length theory. If we restrict consideration to 3D 
compressible convection simulations for simplicity and 
homogeneity, there arc several previous studies which 
have confronted mixing length theory to which we can 
compare our results. These studies investigate convec- 
tion under diverse conditions, including slab convection 
(Chan & Sofia 1987, 1989, 1996; Porter & Woodward 
2000), a red giant envelope (Porter, Woodward, & Ja- 
cobs 2000), and solar and sub- giant surface layers (Kim 
et al. 1995, 1996; Robinson et al. 2003, 2004, 2005). The 
number of zones used range from 1.9 x 10^ (Chan & Sofia 
1989) to 6.7 X 10'^ (Porter & Woodward 2000). The equa- 
tions of state used include a gamma-law (Chan & Sofia 
1989; Porter & Woodward 2000), ionized gas (Kim et 
al. 1996; Robinson et al. 2004), and a combined rela- 
tivistic electron plus ion gas (Timmes k. Swesty 2000) 
in this paper. Subgrid scale physics was treated by a 
Smagorinsky model (Smagorinsky 1963) or by ignoring 
it. We note that Styne, et al. (2000) have shown that 
PPM methods solving the Euler equations converge to 
the same limit as solutions to the Navier-Stokes equa- 
tions, as resolution is increased and viscosity reduced. 
Additionally, the subgrid scale turbulence "model" im- 
plicit in the numerical algorithm of PPM is known to 
be well behaved (Fernando Grinstein, personal commu- 
nication). Given this already inhomogeneous set of sim- 
ulations, determining consistent convection parameters 
is difficult. Our attempt is given in Table 4, in which 
we summarize the convection parameters found in these 
studies for comparison to our own. 

How well do these compare? In some respects the 
agreement is striking. The parameter as is in the range 
^0.7 - 0.8 for all groups. Further, all agree that for their 
case, the mixing-length theory gives a fairly reasonable 
representation of the simulations in the sense that the 



alphas are roughly constant throughout the body of the 
convection zone. The difficulty is that the specific val- 
ues of these alphas depend upon the case considered. 
The two best-resolved simulations, ours and (Porter & 
Woodward 2000), use the same solution method, PPM, 
yet have the most differing alphas. This suggests to us 
that the differences arc due to the physical parameters 
of the respective convection zones. Porter, Woodward, & 
Jacobs (2000) have already shown that slab and spheri- 
cal geometry give qualitatively different behavior for the 
alphas. Our shell is only two pressure scale heights in 
depth, and is relatively slab-like; Porter & Woodward 
(2000) have a convection zone which is more than twice 
as deep by this measure. There is a suggestion in Table 4 
that the alphas increase with the depth of the convection 
zone. This would be reasonable if a convective plume 
were accelerated throiigh the whole convection region 
before it is decelerated at the nonconvective boundary. 
However, the other differences mentioned above proba- 
bly contribute to the scatter in the alpha parameters in 
4. 

Further efforts on this issue are needed. If convec- 
tion does depend upon the nonlocal, physical structure 
of the star, calibration of the mixing length to fit the sun, 
as is traditionally done, is not wise. Furthermore, it is 
well known that the mixing length theory is particularly 
prone to problems in the surface layers where convection 
becomes inefficient. Therefore, the empirical agreement 
of mixing length calibration to the sun and to Popula- 
tion II giants (Ferraro, et al. 2006) may be a fortuitous 
coincidence. 

7. MIXING AT CONVECTIVE BOUNDARIES 

The boundaries which separate the convective regions 
from the stably stratified layers in our 3D simulations 
span a range of relative stability, with 1 < Ris ^ 420. 
At the lowest values oi Ris, the boundary is quickly over- 
whelmed by turbulence, as described in §4.1. Once Ris 
becomes large enough, the boundary becomes stabilizing 
and evolves over a much longer timescale. Snapshots of 
the quasi-steady shell burning and the core convection 
boundaries are presented in Figure 23, ordered by Ris 
with spans the range 36 < Ris ^ 420. The anatomy 
of the convective interfaces includes the turbulent con- 
vection zone, the distorted boundary layer of thickness 
h, and the stably stratified layer with internal wave mo- 
tions (compare to Figure 1). 

The boundary becomes more resilient to thickening, 
and distortion by the turbulence as Ris increases. A 
region of partial mixing exists primarily on the turbu- 
lent side of the interface, where material is being drawn 
into the convection zone. The "ballistic" picture of pen- 
etrative overshooting (Zahn 1991) in which convective 
eddies are envisioned to pierce the stable layer does not 
obtain. Instead, material mixing proceeds through in- 
stabilities at the interface, including shear instabilities 
and "wave breaking" events, which break the boundary 
up into wisps of material that are then drawn into the 
turbulent region and mixed. The convective interface re- 
mains fairly sharp in all cases, and the effective width 
is well described by the elastic response of the bound- 
ary layer to incoming eddies, h ^ Vc/N. The convective 
interfaces seen in our simulations bear a striking resem- 
blance to those observed in laboratory studies of turbu- 
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lent entrainment of comparable Ris (see e.g. McGrath 
et al. 1997, Figs. 2-5). 

The mixing that occurs due to the instabilities and 
eddy scouring events at the interface leads to a steady 
increase in the size of the convection zone. In this sec- 
tion we quantify the entrainment rates at the convective 
boundaries, wc discuss these results in terms of the the 
buoyancy evolution of the interface, and wc describe how 
the "turbulent entrainment" process can be incorporated 
into a stellar evolution code as a dynamic boundary con- 
dition to be used in addition to the traditional static 
Lcdoux and Schwarzschild criteria. Wc conclude the sec- 
tion with some comments on numerical resolution. 

7.1. Quantifying the Boundary Layer Mixing Rates 

As evident in Figure 23, the convective boundary lay- 
ers arc significantly distorted from spherical shells. To 
estimate the radial location of the interface we first map 
out its shape in angle ~ ri{9,(j)). At each angular posi- 
tion the surface is taken to be coincident with the radial 
position where the composition gradient is the steepest 
(this is comparable to the location of minimum density 
scale height Hp = [dlnp/dr]~^). The interface thickness 
h is taken to be the r.m.s. variation of the surface rj with 
angle, h = cr[ri{9, (/>)], which provides a quantitative mea- 
sure of the amplitudes of the distortions imparted to the 
interface. The mass interior to the interface is calculated 
according to, 

{r^) 

= y" 4Trr^{p)dr (27) 

ro 

where ro is the inner boundary of the computational do- 
main, (p) is the horizontally averaged density, and the 
mean interface radius is used for the upper limit on the 
integral. The time derivative Mi is the rate at which 
mass is entrained into the convection zone. 

In Figure 24 the time histories of the averaged interface 
location (r,) and interfacial thickness h are shown for the 
convective boundaries in our simulations. A 3D model 
and a representative 2D model are shown for each bound- 
ary. The outer shell boundary layer adjusts rapidly in the 
first 100s to a new position, due to the penetration event 
discussed in §4.1, after which a slow outward migration 
ensues. For the 3D model, the outward migration pro- 
ceeds in distinct stages, labeled (a - c) in Figure 24. Each 
stage is well described by a linear increase of radius with 
time, and ends with a rapid adjustment to a new entrain- 
ment rate. This behavior can also be seen in Figure 4, 
where the change in entrainment rate is coincident with 
changes in the background composition gradient and sta- 
bility (compare to the initial buoyancy frequency profile 
in Figure 2). 

The downward migration of the lower shell boundary 
is more uniform and proceeds at a significantly reduced 

rate compared to the upper boundary. The core convec- 
tion boundary evolution departs most significantly from 
a linear trend, but Monotonic growth is clearly estab- 
lished very soon after the simulation begins t>2x lO^s. 

The interfacial thickness h in the oxygen burning mod- 
els are initially large, reflecting the strong mixing event 
during the initial transient, but settle down to relatively 
constant values for t > 300s. In contrast, the boundary 



thickness in the core convection model increases gradu- 
ally with time until a steady state value is achieved, re- 
flecting the milder initial development. In all cases, the 
time averaged values of h during the quasi-steady states 
compare well to the boundary displacement expected for 
eddies impacting the stable layer with the characteristic 
convective velocities of the model, h ^ Vc/N. 

The entrainment rate and the interfacial thickness is 
larger in all of the 2D models as a consequence of the 
larger velocity scales in those simulations. The interface 
migration rates and averaged interfacial layer thicknesses 
are tabulated for all of the models in Tables 5 and 6, and 
are broken down into various time intervals over which 
linear growth of the boundary is a good approximation. 
Time averaged mass entrainment rates are also included 
in Tables 2 and 3. 

7.2. The Entrainment Energetics 

In order for entrainment to take place at a convec- 
tive boundary the buoyancy increment of the stable layer 
material over that of the mixed layer material must be 
overcome. This can happen in two distinct ways. First, 
non-adiabatic processes can change the relative stability 
of the stable layer. For example, heating the convective 
region will cause an increase in its entropy, and the buoy- 
ancy jump separating the overlying layer will decrease. 
The rate at which the convective boundary will grow due 
to heating is Us = s/{drs), where drS is the radial gra- 
dient of entropy at the boundary and s is the time rate 
of change of entropy in the shell. This process will cause 
both the upper and lower boundaries to migrate to larger 
radii the upper boundary will be weakened, while the 
lower boundary will become stiffer. Non-adiabatic pro- 
cesses in the boundary layers will affect their stability in 
the same way: cooling in the upper and heating in the 
lower boundaries will weaken their stratification. 

A related, but distinct process is "turbulent en- 
trainment" whereby turbulent kinetic energy does work 
against gravity to draw material into the turbulent re- 
gion. In this process, the stratification is weakened at 
a convective boundary by the turbulent velocity fluctu- 
ations. This is quantified in terms of the buoyancy fiux 
q = gp'v' /po- In the absence of heating and cooling 
sources the buoyancy in the interfacial layers will evolve 
according to the buoyancy conservation equation, 

dtb = -div{q) (28) 

and a positive flux divergence at the boundary will lead 
to a weakening of the stratification. The relationship 
between turbulent entrainment and the weakening of a 
boundary through heating and cooling can be understood 
in terms of the enthalpy flux which attends the buoyancy 
flux. In fact, the buoyancy fiux is directly related to the 
enthalpy flux across the interface, 

= poCp(TX) = ^(p'O = pocp^ X q (29) 

and is equivalent to heating and cooling processes oper- 
ating in the boundary layer (note the downward directed 
enthalpy flux within the boundary layers in Figure 15). 

What drives the entrainment seen in the present sim- 
ulations? Can the entrainment in the outer shell bound- 
ary be explained by the heating of the convection zone 
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by nuclear burning? Comparing the entropy growth rate 
of the shell to the entropy gradient at the boundary we 
find Us ~ 0.8 x 10'' cm/s, which is at most 17% of the 
measured growth rate for this boundary, and typically 
of order a few percent. Shell heating will reduce the 
growth rate of the lower boundary by ~ 0.04 x 10^ 
cm/s, which is of order a few percent of the rate mea- 
sured. Therefore, the overall heating and cooling of the 
shell contributes very modestly to the growth of the shell 
over the course of the simulation. The long thermal 
timescale in the core convection model reduces this effect 
even more, where it is lower by several orders of magni- 
tude. Therefore, we turn to the turbulent hydrodynamic 
processes operating in the boundary layer to understand 
the growth of the convection zones. 

In Figure 25 we present the buoyancy flux profiles for 
our 3D simulation models, including both time-series di- 
agrams and time averaged radial profiles. The properties 
of the buoyancy flux can be divided into three distinct 
flow regimes: (1) the body of the buoyant convecting 
layer, which is dominated by positive q; (2) the con- 
vective boundary layers, with negative q; (3) the stably 
stratified layers, where q is oscillatory, but has a nearly 
zero mean (in both a horizontal and time average sense). 

The buoyancy driving of the convective flow in regime 
(1) can be appreciated by comparing the flow velocity 
to the commonly used buoyant convection velocity scale 
■= 2.5 J {q)dr, where integration is taken over the ra- 
dial extent of the convection zone (see e.g. Deardorff 
1980). for the 3D shell burning and core convection 
models are ~ lO^cm/s and t;, ~ 3 x 10^ cm/s, which 
compare well to the radial r.m.s. turbulent velocities 
measured in the simulation (Figures 6 and 12). 

In regime (2), which occurs in the convective bound- 
aries, the buoyancy flux is negative. A negative value of 
q signifies that the turbulent kinetic energy is being con- 
verted into potential energy. The mixing that attends 
this negative buoyancy flux underlies the entrainment 
that is taking place at the boundaries through equation 
28. We demonstrate this by showing that the entrain- 
ment speeds measured in the simulation are consistent 
with the measured buoyancy fluxes. The interface migra- 
tion speed is incorporated into the conservation equation 
by writing the time derivative as an advective derivative, 

dtb ~ Uedrb = UeN'^ (30) 

where we have used the relationship drb = N"^. Using 
this time derivative in equation 28 and solving for Ug we 

find. 



where we have approximated the divergence of the buoy- 
ancy flux with the difference Aq/h. We use the symbol 
Ug to distinguish the estimated rate from the values mea- 
sured in the simulation. 

If we adopt the buoyancy flux at the interface for Aq 
(Figure 25), the measured interface thickness h, and the 
buoyancy frequency at the boundary, we flnd the follow- 
ing entrainment rates. For the upper shell boundary, 
lower shell boundary, and the core convection boundary 
we have: Ue ~ 5.1 x lO'' cm/s; Ug ~ 1.1 x 10^ cm/s; 
and Ue ^ 2.2 x 10^ cm/s, respectively. These are to be 



compared with Ug — jr^ — Vg^pl measured in section §7.1 
and presented in Tables 5 and 6. The values correspond- 
ing to the same time period are: Ug = 4.6 x 10** cm/s; 
Ug = 1.2 X 10^ cm/s; and Ug = 2 x 10^ cm/s. Although 
these estimates are only order of magnitude (e.g., using 
the crude approximation for the time derivative in eq. 
[30]) they compare well to the values measured in the 
simulations and the buoyancy flux due to " turbulent en- 
trainment" can account for the growth of the convective 
layers seen here. 

7.3. Whence q? 

The buoyancy flux q appears as a term in the turbu- 
lent kinetic energy (TKE) equation, which we present in 
§A (eq. [A12]). In our notation, the buoyancy flux is re- 
lated to the buoyancy work term by g = {Wb)/po- The 
buoyancy flux, therefore, is related to the rate at which 
turbulent kinetic energy is advected into the stable layer 
Fk, the rate at which it dissipates through viscous forces 
Ek, and the rate at which energy is transported through 
the boundary layer by pressure- velocity correlations Fp. 
In essence, entrainment is the process by which the tur- 
bulent kinetic energy in the boundary layer does work 
against gravity to increase the potential energy of the 
overall stratification. 

Two theoretical approaches have been taken to study 
entrainment. The first approach ignores the TKE equa- 
tion and instead posits an " entrainment law" . The en- 
trainment law is merely a functional form for the rate 
at which stable layer mass will flow into the turbulent 
region, and is therefore a dynamic boundary condition. 
These laws are visually parameterized by the stability 
properties of the interface and the strength of the tur- 
bulence through RiB (see e.g. Fedorovich, Conzemius & 
Mironov 2004). Once an entrainment law is adopted, 
the enthalpy flux can be calculated and the evolution of 
the boundary can be self-consistently solved for. The 
advantage of such an entrainment law is the simplicity 
with which it can be incorporated into global circulation 
models of the atmosphere, for instance. 

An alternative approach to adopting an entrainment 
law is an explicit physical model for the terms in the TKE 
equation (eq. [A12]). For example, general forms for the 
buoyancy flux profile within the stable layer have been 
applied with some success in reproducing the growth of 
the atmospheric boundary layer and the deepening of 
the oceanic thermocline (StuU 1976b; Deardorff 1979; 
Fedorovich & Mironov 1995). In some respects, how- 
ever, these models are glorified entrainment laws since 
the buoyancy flux is prescribed in a simplifled, parame- 
terized way. Moving beyond assumptions concerning the 
turbulence profiles within the interfacial layer, are theo- 
retical models which take into account the interactions of 
waves and turbulence and incorporate non-linear models 
for the evolution of instabilities (e.g. Carruthers & Hunt 
1986; Fernando & Hunt 1997). The approach adopted 
in these theoretical studies is general enough that any 
adjustable parameters may turn out to be universal and 
a prcKlictivc model can be developed. In addition, the 
framework employed is general enough that the produc- 
tion of turbulence by mean flows (i.e., stellar rotation) 
can be incorporated, as well as long-range effects due 
to internal waves. The internal waves are incorporated 
through the pressure-correlation flux, Fp, and plays a 
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central role in the evolution of the buoyancy flux when 
wave breaking is important. 

7.4. An "Empirical" Entrainment Law 

The development of a sophisticated turbulence model 
to explain entrainment is beyond the scope of the present 
work. Instead, we ask to what extent do the entrainment 
laws used in geophysical models apply to our simulations 
and stellar interiors? Guided by laboratory study and 
geophysical large eddy simulation we study the depen- 
dance of the entrainment rate on the bulk Richardson 
number. 

Rig is calculated according to equation 1, where wc use 
the horizontal correlation length scale L = Ch defined in 
§B. The buoyancy jump is calculated by performing the 
integration in equation 2 over the width of the interface 
which we take to be the interval r G [Ti — h,ri + h]. The 
normalized entrainment rates E = Ue/cr, the buoyancy 
jumps Ab, and Ris are presented in Tables 5 and 6. The 
dependance of the entrainment coefficient E on Ris is 
presented in Figure 26. 

The 2D and 3D data are found to obey similar trends 
(lower E for higher Rib), but occupy signifanctly differ- 
ent regions of the diagram. This can be explained by the 
much higher r.m.s. velocities in the 2D simulation. The 
velocity scale in 2D is apparently an artifact of the re- 
duced dimensionality of the problem which significantly 
influences the flow morphology. Although the velocity 
scale is higher in the 2D models, it is much more laminar 
and accompanied by less turbulent mixing. The arrow in 
Figure 26 indicates the direction that the 2D data points 
would move if a lower effective r.m.s. velocity were as- 
sumed. In what follows we focus our attention exclusively 
on the entrainment data found for the more realistic 3D 
models. 

What we find is that the entrainment coefficient E is 
well described by a power law dependance on Ri b of the 
form in equation 4. Our best fit values for the parame- 
ters are log A = 0.027 ± 0.38 and n = 1.05 ± 0.21. This 
entrainment law is shown by a dashed line in Figure 26. 
Remarkably, the power law is of order unity, in agree- 
ment with geophysical and laboratory studies. The fact 
that the entrainment in our simulations are governed by 
the same, fairly universal dependance on RiB as these 
other studies may have been anticipated, considering the 
striking degree of similarity between the buoyancy pro- 
files and the character of the developed fiow in the vicin- 
ity of the boundary (Figure 23). 

7.5. A Dynamic Convection Zone Boundary Condition 

Mass entrainment is a fundamentally different phe- 
nomena from diffusion, which is the typical route used 
to incorporate new mixing phenomena into a stellar evo- 
lution code. Therefore, how might we incorporate this 
new process? Schematically, the idea is very simple. For 
each convcctive boundary, initially found with the tradi- 
tional stability criteria {ds/dr = 0, d^s/dr^ 7^ 0), we can 
calculate the associatcid bulk Richardson number based 
on the background stratification and an approximation 
of the turbulence characteristics (e.g., from MLT). With 
Ri B in hand we can then input this into our entrainment 
law, E = E{RiB) which returns to us the entrainment 
rate. The entrainment rate, therefore, is the boundary 



growth rate as a function of Ri b and possibly other pa- 
rameters of the system. The function ElRis) can be 
broken up into at least three regimes for convenience. 

Low stability: RiB < Ri'g^"'. For low RiB it is ob- 
served that mass entrainment happens very quickly, on 
an advection timescale (§4.1). Therefore, we can define 
a minimum _R«^*" at which the expansion of the con- 
vection zone will proceed very quickly, eliminating con- 
vcctive boundaries which are too weak to support the 
adjacent turbulence. 

Intermediate Stability: Ri'g''' < RiB < Ri']^"''. For 
an intermediate range of stability, we can use the fairly 
universal entrainment law which matches our simulation 
data, defined by the two parameters A and n. Although 
scatter in mixing rates were found to be as large as a 
factor of a few relative to the best fit law, the general 
monotonic, power-law dependance was found to be ro- 
bust. We can incorporate this physics into the stellar 
evolution code as a mass entrainment rate. 

Me = = {^nr^Pi)cTH x /a x loi-niosRis) (32) 

where the normalization factor is written = 10*^'°^"*^ 
and represents the turbulent entrainment mixing effi- 
ciency. More sophistication can subsequently be incor- 
porated as our understanding of the entrainment process 
improves. 

High Stability: RiB > Ri^""^. The entrainment pro- 
cess will cease to operate at some upper limit Ri^"''^ , 
above which the boundary evolution will be controlled 
by diffusive processes on the molecular scale. Following 
(Phillips 1966), we have, 

«^™=(^)(^) (33) 

which is based on the condition that the kinetic energy 

in the turbulence is sufficient to lift the material from 
the interface, pa^ > pN'^A?. Here, the interface thick- 
ness is taken to be that due to molecular diffusion with 
A > k/ue- The relatively small diffusion rates in stel- 
lar interiors imply that turbulent entrainment will con- 
tinue to operate to very high Richardson numbers. For 
comparison, the entrainment process in the ocean is es- 
timated to operate up to Rib ^ 10^^®. 

Additional details concerning the implementation of 
this type of boundary condition into TYCHO will be 
presented in a subsequent paper. 

7.6. Spatial Scales, Numerical Resolution, and 

Entrainment 

We conclude this section with a few comments on how 
well simulation can be trusted in elucidating the process 
of entrainment, which is not very well understood. The 
spatial scales which limit the entrainment rate at a con- 
vcctive boundary are not also not well understood, and 
depend on the interplay between large eddy and small 
scale turbulent transport processes. As discussed by 
Lewellen & Lewellen (1998), there is feedback between 
the transport rate away from the turbulent boundary 
layer which is controlled by large eddies, and the trans- 
port rate of material in the immediate vicinity of the 
interface by small scale turbulence. A full understanding 
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of this problem hinges on being able to resolve the en- 
trainment zone in the presence of the large scale eddies. 

A code comparison and resolution study of the entrain- 
ment problem in the planetary boundary layer context 
was conducted by Bretherton et al. (1999) and Stevens 
& Bretherton (1999). In these studies, it was suggested 
that the appropriate criteria for resolving boundary layer 
entrainment is that the grid zoning is smaller than the 
fluctuations induced in the inversion layer by the large 
eddies so that shear instabilities (e.g., Kelvin-Helmholtz) 
would not be suppressed. It was suggested that a "nested 
grid" of refinement within the boundary layer was com- 
parable to using fine resolution throughout the simula- 
tion domain. The suggested resolution criterion in these 
papers, however, fail to account for the non-linear turbu- 
lent evolution which proceeds the onset of instabilities. 
In addition, no simulations were conducted with enough 
resolution that the boundary layers were turbulent, and 
the simulations presented were even marginally resolved 
by the author's suggested criteria. 

A related study by Alexakis (2004) investigates the en- 
trainment and mixing at a boundary due to internal grav- 
ity wave breaking driven by a shear flow. This process 
may be operating in the shear mixing layers that form 
when large eddies impact the boundary layer. In this 
study, the mass entrainment rate was found to depend 
on the numerical resolution in a non-monotonic way, first 
decreasing and then increasing with finer resolution. The 
author concludes that low resolution models are domi- 
nated by numerical diffusion until the resolution is fine 
enough to resolve turbulence near the boundary, at which 
point the entrainment rate begins to increase and is con- 
trolled by turbulent transport. Although the asymptotic 
mixing rate was inconclusive and no resolution criteria 
was proposed, resolving the turbulence ensuing from the 
instability was shown to be important. 

Much more work needs to be done to address the role 
played by both the small scale processes and the inter- 
play with large eddies. Two complimentary numerical 
approaches can be taken. First, the feedback between 
large and small scale mixing and transport processes 
can be studied using large eddy simulation with a range 
of subgrid scale mixing efficiencies. Such a study can 
help develop insight into which scales control the mixing 
rate. Second, direct numerical simulations (DNS) which 
resolve the turbulent processes operating at the inter- 
face can be undertaken when sufficient computational re- 
sources are available. These studies would provide more 
definitive conclusions concerning the interplay between 
eddy scales and would provide guidance for a more gen- 
eral framework for future theoretical analysis. Finally, it 
is important to keep in mind that laboratory studies of 
high Reynolds number turbulent entrainment continue 
to provide useful constraints, and improved flow visual- 
ization techniques are allowing a more direct comparison 
to theory and simulation. 

The "empirical" entrainment law which we discuss in 
this paper is constrained by only a few data points (the 
six 3D data points in Figure 26). Extending simiflations 
to include an ever more diverse suite of stellar structures 
would provide an even stronger mandate, and better con- 
strained model for incorporating this physics into stellar 
evolution codes. 



8. SUMMARY AND CONCLUSIONS 

In this paper, we have presented the results of three- 
dimensional, reactive, compressible, hydrodynamic sim- 
ulations of deep, efficient stellar convection zones in mas- 
sive stars. Our models are unique in terms of the degree 
to which non-idealized physics have been used, and the 
evolutionary stages simulated, with fuel and ash clearly 
distinguished. 

We find several general results regarding the basic 
properties of the convective flow: 

• the flow is highly intermittent, but has robust sta- 
tistical properties, 

• the 2D vs 3D velocity scales differ by almost a fac- 
tor of several, and the flow morphologies are com- 
pletely different, 

• stable layers interact with convection to deceler- 
ate plumes, and consequently distort these layers, 

which then generate waves, 

• mixing is found to occur at convective boundaries 
in manner best described as turbulent entrainment, 
rather than the traditional picture of convective 
overshooting wherein turbulent eddies ballistically 
penetrate the stable layers. 

We have compared our oxygen shell burning model to 
mixing length theory assumptions. We show that, while 
a reasonable representation of the super-adiabatic tem- 
perature gradient and velocity scale can be fit with a 
single mixing length, the values of the inferred mixing- 
length "constants" differ from other simulations. This 
was already implied in Porter, Woodward, & Jacobs 
(2000), who found difference for slab and spherical ge- 
ometries. There may be a dependence upon the depth of 
the convection zone as well, and possibly upon the na- 
ture of the stable boundary regions and/or the nature of 
the driving process (burning or radiative loss). 

Why do we care about MLT in regions of efficient con- 
vection? (1) temperature profile can affect the burning 
rates, which have a stiff temperature dependence; (2) 
the velocity scale can affect the nucleosynthesis (such as 
s-process branching ratios in double shell burning AGB 
stars) by dictating the exposure time of the plasma to 
varying conditions throughout the burning region; (3) 
the velocity scale and the kinetic energy flux is an im- 
portant input needed for calculating the mixing at con- 
vective boundaries. 

We have found that the extent of mixing is better 
represented by an integrated Richardson number rather 
than the convectional Schwarzschild or Ledoux criteria 
alone. This incorporates the addition physics related to 
the resistance of stiff boundaries to mixing. Related to 
the definition of boundary stiffness, we have identified 
an important physical process which is missing from the 
standard theory of stellar evolution: turbulent entrain- 
ment. This process is well known in the meteorology and 
oceanographic communities, and has been extensively 
studied experimentally. We show that the rate of en- 
trainment is well represented as a simple function of the 
buoyancy jump, in a manner similar to that determined 
experimentally. 
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The long term consequences of convectivc boundary 
inconsistencies such as the one illustrated by the initial 
transient in our simulation, and for which the conditions 
are common in ID stellar models, can significantly al- 
ter the size of convective cores, and thus the subsequent 
explosion and nucleosynthetic yields of the resultant su- 
pernova. In a subsequent paper in this series, we will 
present case studies which incorporate the physical in- 
sight gained through these simulations into the TYCHO 
stellar evolution code. We expect to see effects in solar 
models, s-processing in AGB stars, stellar core formation 



(white dwarfs, neutron stars, and black holes), stellar nu- 
cleosynthesis yields, stellar ages, and HR diagrams. 
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APPENDIX 
THE ENERGY EQUATION 
Total Energy 

The primitive energy equation solved by PROMPI is. 



dt {pE) + V • [{pE + p)u + Fr] = pu • g + pCnet 



(Al) 



where the total energy is composed of the internal and kinetic components, E = Ei + Ek- We decompose the velocity, 
density, and pressure fields into mean and fluctuating components according to, 



where {ip) = <po and {ip') = 0. The pressure-velocity correlation term is. 



The gravity term is, 



V • (pu) = V • (pouo) + V • (poW) + V • (p'uo) + V • (p'u'). 



(pg • u) = (pouog) + (pou'g) -I- (p'uog) -I- (p'u'g). 
The averaging operator eliminates terms which are first order in fluctuations (by definition) and we have. 



dtipE) 



(pEuo) + {pEu') + (pouo) + {p'u') 



(Pouog) + (p'n'g) + (penet)- 



(A2) 



(A3) 



(A4) 



(A5) 



We can further simplify this expression using the condition of hydrostatic equilibrium, which holds to a high degree of 
accuracy in the simulation (Vpo = Pog)- The background velocity in this case, Uq, is a slow, highly subsonic, expansion 
or contraction that is driven on a thermal timescale. The background velocity field has only a radial component (i.e., 
there is no rotation in the current model), Uq = (uo,(j.), 0, 0). The energy equation can be then simplified to read, 



dtipE) + V ■ (pEuo) =-V 
where we have used the following definitions, 



(Fp + Fi + Fk + Fr) - (poV ■ Uo) + (Wb) + {penet) ■ 



(A6) 



Fi 

Fk = 



= pEiu' 

pExu' 



bp =pu 
Wb = p'g ■ u'. 



(A7) 
(A8) 
(A9) 
(AID) 



Kinetic Energy 

The kinetic energy equation is derived by forming the scalar product of the velocity with the equation of motion 
(e.g., Shu, 1992, Ch.2). The momentum equation can be written in vector form as, 



dtipEx) + V • {pExu) = -u • Vp -I- pu • g 



(All) 



Again, we decompose the fields into mean and fluctuating components, employ the hydrostatic equilibrium condition, 
and perform averages. The result is. 



dtipEK) + V • {pEkUo) = -V • (Fp + Fk) + (p'V • u') + (Wb) - sk- 



(A12) 



Here, Ek is the viscous dissipation of kinetic energy. In our simulations, this term is not modeled explicitly and arises 
due to numerical dissipation. The term p'V • u' represents the compressional work done by turbulent fluctuations, and 
the other terms are as deflned above. 
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CORRELATION LENGTH SCALES 

The vertical correlation of the horizontal distribution of fluctuations in a quantity X' = X — {X) a.t radial position 
r and offset position r + 6r is calculated according to, 

^ ^ An (Tx{r)crx{r + Sr) ^ ' 

where the integral is taken over the angular direction with dfi = sm{9)d9d<j). The correlation is normalized by the 
product of the horizontal r.m.s. value of the quantity at the two levels being compared ax- 
The horizontal correlation of fluctuations at radial position r is calculated using the autocorrelation function, 

where the brackets (•) denote averaging over all horizontal locations s and fixed offset Ss. The horizontal correlation 
is normalized by the variance of the quantity cr^ . 

Characteristic length scales are defined as the offset position where the correlation function drops to a value of 0.5. 
For horizontal correlations, we define this length as jCh- We also define a value which is twice this length, the full 
width at half maximum, which we denote by Lu- (The value Ch provides a good approximation to the integral scale, 
JC"{6s;r)d6s.) 

In the vertical direction the sign of the offset 6r is retained and a separate length scale is defined where the correlation 
function drops to 0.5 for positive and negative offsets, which we denote by Ly and Ly. The full width is denoted 

Lv = Ly — Ly. 
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Fig. 1. — Diagram illustrating the salient features of the density and velocity field for the turbulent entrainment problem. Three layers 
are present: a turbulent convection zone is separated from an overlying stably stratified region by a boundary layer of thickness h and 
buoyancy jump Ab ~ N'^h. The turblence near the interface is characterized by integral scale and RMS velocity Ch and cr/f, respectively. 
The stably stratified layer with buoyancy frequency N{r) propagates internal waves which are excited by the adjacent turbulence. A shear 
velocity field v_i_{r), associated with differential rotation, may also be present. After Strang & Fernando (2001). 
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Fig. 2. — Radial profile of the simulated region for the oxygen shell burning models. The thin lines indicate the initial conditions and 
the thick lines indicate the 3D model at t = 400 s. (top left) Temperature and density, (top right) Mass fraction of ^®0. (bottom left) 
Squared buoyancy frequency, (bottom right) Buoyancy. 
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Fig. 3. — This time sequence shows the onset of convection in the oxygen shell burning model. The first 200 s of the 2D model (ob.2d.c) 
is shown, including the initial transient and the settling down to a new quasi-steady state. The light yellow line indicates the location of 
the convective boundary as defined in the ID TYCHO stellar evolution model (Ledoux criterion), which was used as initial conditions for 
the simulation. 
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Fig. 4. — The time evolution of the 3D oxygen shell burning model, (top) The magnitude of the oxygen abundance gradient is shown 
and illustrates the migration of the convective boundaries into the surrounding stable layers. Interfacial oscillations are also apparent in 
the upper convective boundary layer (r ~ 0.85 X 10*' cm), and internal wave motions can be seen quite clearly in the upper stable layer, 
(bottom) The kinetic energy density is shown, and illustrates the intermittent nature of the convective motions. The upwelling chimney-like 
features in the convective region are seen to excite internal wave trains in the stable layers, which propagate away from the boundaries of 
the convection zones. See also Fig. 25. 
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Fig. 5. — The time evolution of the energy budgets for the oxygen shell burning models: the (thick line) 3D model, and (thin lines) the 
three 2D models are shown, including: (thin-solid) ob.2d.c; (thin-dashed) ob.2d.e; and (thin dotted) ob.2d.C. The energy budget includes 
the internal energy Ej, the gravitational energy Eq, and the kinetic energy Ex- Note that the energy scale is logarithmic, so that the 3D 
kinetic energy is much smaller than the 2D values. 





Fig. 6. — The r.m.s. velocity fluctuations for oxygen shell burning: (left) 3D model, with velocity components (thick-solid) Vr, (thin- 
solid) vg, and (thin-dashed) v^. (right) The 2D models, with velocity components (thick) Vr and (thin) for simulations (solid) ob.2d.e, 
(dashed) ob.2d.c, and (dash-dot) ob.2d.C. 




Fig. 7. — Mode diagrams for several radial positions in the oxygen shell burning model show the dominant spatial and time scales on 
which motions occur. The abscissa measures k which is related to the wavenumber index / of the mode hy I = 12 X k. The three locations 
shown here include: (left) Lower stable layer, just beneath the convective shell r = 0.4 X 10^ cm. (middle) Middle of convective shell, 
r = 0.6 X 10^. (right) Upper stable layer, just above the convective shell r = 0.9 X 10® cm. 
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Fig. 8. — Radial profile of the simulated region for the main sequence core convection model. The thin lines show the initial conditions 
and the thick lines show the state of the 3D model at t = 10® s. (top left) Temperature and density, (top right) Hydrogen abundance, 
(bottom left) Squared buoyancy frequency, (bottom right) Buoyancy. 




Fig. 9. — Velocity isocontours show the development of the flow in the 3D core convection model. The turbulent convective flow excites 
internal waves which radiate into the overlying stably stratified layer. By the end of the time sequence shown the stable layer cavity is 
filled with resonant modes. 
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Fig. 10. — The velocity magnitude for the core convection model at t=10® s: (left) a slice through the 3D model; and (right) the 2D 
model. The topology of the convective flow is significantly different between 2D and 3D models: the 3D convective flow is dominated by 
small plumes and eddies while the 2D flow is much more laminar, and dominated by a large vortical eddies which span the depth of the 
layer. The wave motions in the stable layer have similar morphology in 2D and 3D, but the velocity amplitudes are much larger in 2D. 
The computational domains have been tiled once in angle for presentation. 




Fig. 11. — The time evolution of the energy budget for the main sequence core convection models: the (thick line) 3D model; and (thin 
line) the 2D model are shown. The energy budget includes the total internal energy Ej, gravitational energy Eq, and kinetic energy Ej^ 
on the computational grid. 




Fig. 12. — The r.m.s. velocity fluctuations for the core convection model: (left) the 3D model, and (right) the 2D model. In each plot, 
the thick line indicates the radial velocity component and thin line is used to indicate horizontal velocity components, with the dashed line 
used to show the polar angle component in the 3D model. 



24 



Meakin & Arnett 




1.0 1.5 2.0 2.5 

Rodius (10" cm) 



Fig. 13. — Pressure fluctuations in core convection model: The time averaged horizontal r.m.s. pressure fluctuations are shown as the 
thick Hne, with extreme values over two convectivc turnovers indicated by the shaded region. The thin line shows the radial component of 
the turbulent ram pressure p'o^ averaged over a convectivc turnover. At the upper boundary, the curves cross at a point where the turbulent 
pressure is balanced by the wave induced pressure fluctuations in the stable layer. This crossing point is coincident with the location of the 
convective boundary. The pressure perturbations at the lower boundary are due to the input luminosity which drives the convective flow. 




Fig. 14. — (left) Density fluctuations in the 3D core convection model; The time averaged maximum density fluctuation is shown as the 
thick line, with extreme values for the averaging period (two convectivc turnovers) shown by the shaded region. The largest fluctuations 
occur at the interface between the turbulent convectivc region and the stably stratified layer. The maximum fluctuation at the interface 
is p'/(p) ~0.12%. (right) The buoyancy frequency is shown in units of (10^* rad/scc). Also shown by the dashed line is the buoyancy 
frequency normalized by the gravity which sets the scale of the density fluctuations at the convective boundary through equation 13. The 
expected density fluctuation is p' /{p) ~ 'Vc\N\/g ~ 0.12%, where a velocity scale of iic ~ 2 x 10® cm/s has been used (see Figure 12). 
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Fig. 15. — (left) Convective enthalpy flux, Fc = {pCpVrT'). (right) Temperature-velocity correlation function as calculated according to 
equation 21, with mean value {cxe) = 0.7 shown by the dashed line. 
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Fig. 16. — (left) Dimensionless temperature gradients: the stellar interior Vs; adiabatic Vo^; a^nd Ledoux Vjed gradients are shown, 
(right) Super-adiabatic temperature gradient horizontally and time averaged. 
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Fk;. 17. (left) Time averaged r.m.s. temperature fluctuations: (thick solid) line shows the r.m.s. fluctuations; the (thin solid) and (thin 
dotted) lines show the fluctuations in the upward and downward directed flow components, respectively, (right) The radial dependence of 
the "thermal mixing length" parameters cc-r defined by equation 22 are shown the temperature fluctuations presented in the left panel, 
using the same line types. The mean values, averaged over r G [0.5, 0.75] x lO^cm are shown by the thin dotted lines. 
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Fig. 18. — (left) Radial velocity amplitudes: (thick solid) r.m.s. value; the (thin solid) and (thin dashed) show the mean up- and down- 
flow velocities, respectively, (right) The radial dependance of the "velocity mixing length" parameters defined by equation 23 are shown 
for the velocity amplitudes presented in the left panel, using the same line types. The mean values, averaged over r e [0.5,0.75] x 10^ cm 
are shown by the thin dotted lines. 



0.4 ; E 0.4 




Fig. 19. — The vertical correlation length scales Ly as defined in §B. (top left) Ly for velocity fluctuations, v'^.; (top right) Ly for 
temperature fluctuations, T' . The pressure scale height Hp and density scale height Hp are shown for comparison, (bottom left) Illustration 
of the relationship between eddy shape and the correlation length scales, Ly and Ly. The grey patches represent the shapes of the eddies 

and the Ly^~ values are measured in the radial direction, away from the horizontal line, (bottom right) Correlations lengths Ly scaled to 
pressure and density scale heights, e.g., ay{vr,Hp) = Ly{vr)/Hp 
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Fig. 20. — The horizontal and vertical correlation length scales, Lh (thick line) and Ly (thin line) are shown for temperature (dashed) 
and velocity (solid) fluctuations. 
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Fig. 21. — The fractional area occupied by the upward flowing material fu is shown as a function of radial position. The downward 
flowing area is /d = (1 — /«) and the dashed line at 1/2 indicates up-down symmetry. 




Fig. 22. — (left) Kinetic energy flux: (thick) line shows the value measured in the simulation averaged over two conveetive turnovers; the 
(thin solid) line shows Ffc calculated according to equation 26; the (thin dashed) line shows Fjf calculated according to equation 26 but 
uses c{v)^ in place of (v^), and a correlation constant of c = 5. (right) Asymmetries in radial velocity: the (thick) line show the skewness 
in the velocity field, 7 = {v^)/a^; the (thin solid) and (thin dashed) lines show the correlations x = where the subscripts u and 

d indicate up- and down-flows,respectively. 
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Fig. 23. — Equatorial slices showing the flow in the vicinity of the convective boundaries in the 3D simulations, ordered by relative 
stability: (row a) upper shell convection boundary, Rig ~ 36; (row b) core convection boundary. Rig ~ 48; (row c) lower shell convection 
boundary. Rig ~ 419. The colormap indicates composition abundance, where the darker tones trace stable layer material entrained across 
the interface. The velocity vectors have been sampled every third zone in each dimension. 
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Fig. 24. — The time history of (top row) the convective boundary location, and (bottom row) the thickness of the convective interface 
for: (left) upper shell burning boundary; (middle) lower shell burning boundary; and (right) core convection boundary. The (thick line) 
identifies the 3D models, ob.Sd.B and msc.Sd.B; and the (thin line) identifies the 2D models, ob.2d.e and msc.2d.b. The (dashed lines) 
show the averaged interface thickness for t > 300 s for oxygen burning, and t > 6.0 X 10^ s for core convection. The letters (a-c) in the 
upper left panel mark times when the outward migration rate of the convective boundary rapidly adjusts to a new value in the 3D model. 
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Fig. 25. — Buoyancy flux. Time-series diagrams and time-averaged radial profiles are shown for: (top-row) the 3D oxygen shell burning 
model; and (bottom-row) the 3D core convection model. 



30 Meakin & Arnett 




-4 



10 



100 



Fig. 26. — Normalized entrainment rate plotted against bulk Richardson number Rib- The 3D models are marked with squares, and the 
2D data by plus signs. The best fit power-law to the 3D model data is shown by the dashed line. The 2D entrainment rates fall everywhere 
below the 3D trend. The arrow indicates the direction in the diagram that the 2D data points would move if the effective r.m.s. turbulence 
velocity were lower. 
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TABLE 1 

Nuclei Included in Reduced 

Nuclear Reactiox Network 

Element Charge Atomic 
Weight 

HcHum 2 4 

Carbon 6 12 

Oxygen 8 16 

Neon 10 20 

Sodium 11 23 

Magnesium ... 12 24 

Silicon 14 28 

Phosphorus ... 15 31 

Sulfur 16 32, 34 

Chlorine 17 35 

Argon 18 36, 38 

Potassium .... 19 39 

Calcium 20 40, 42 

Titanium 22 44, 46 

Chromium ... 24 48, 50 

Iron 26 52, 54 

Nickel 28 56 

Note. — Network also includes elec- 
trons, protons, and neutrons. 



TABLE 2 

Summary of Oxygen Shell Burning Models 



Paiauirl rr 


I ; nils 


()l).:id.L5 


()li.2d.c 


()li.2d.(' 


ul>.2<Le 


^ini ^out 


(10« cm) 


0.3, 1.0 


0.3, 1.0 


0.3, 1.0 


0.3, 5.0 


A9,A0 


(deg.) 


30, 30 


90, 


90, 


90, 


Grid Zoning 




400x (100)2 


400x320 


800x640 


800x320 


tmax 


(s) 


800 


574 


450 


2,400 


fconv 


(10'' cm/s) 


0.8 


2.0 


1.8 


1.8 


tconv 




103 


40 


44 


44 


Mi\J 


(lO-^Mo/s) 


1.1 


1.33 


1.25 


1.3 


Mill 


(lO-^Mo/s) 


-0.23 


-0.52 


-0.5 


-0.5 



"'^The subscripts u and / refer to the upper and lower convective shell boundary. 



TABLE 3 

Summary of "Core Convection" Models 

Parameter Units msc.3d.B msc.2d.b 

rf„, lout (10" cm) 0.9, 2.5 0.9, 2.5 

Ae,A4> (deg.) 30, 30 30,0 

Grid Zoning - 400x(100)2 400x100 

tmax (s) 2.0x10*^ 2.0x10^ 

fconv (105 cm/s) 2.5 13 

tconv (s) 3.6 xlO^ exlO"" 

Mi (IO-'^Mq/s) 2.72 4.73 
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TABLE 4 

Assumed and Measured Convection Parameters^ 



Study Pe'' as OiA,T CKA.u ctA L/Hp'^ Grid Zoning 



MLT 

This Studyi 

Chan-Sofia^ 

Kim3 

Robinson* 

Porter- Woodward^'' 



>1 
> 10^ 



10 - 8 X 10* 



0.70±0.03 

0.81±0.03 

0.80±0.01 

0.65-0.85 

0.7-0.9 



a 

0.8 - 
1.32 ■ 
2.96 ■ 

4.08 



1.4 

■ 3.75 

■ 4.2 



a 

1.35 
3.39 
1.5 - 

3.82 



- 1.73 

- 6.4 
3.4 



a 

0.87 
1.90 
1.4 - 



- 1.31 

- 4.4 
3.2 



2.68(3.53)5'' 



~2(3.7)'i 

7 

6 

8.5 

5 



100^ X 223(400)^^ 
20^ X (< 50) 
323 

1142 X 170 
512^ X 256 



^See §6.2 for parameter definitions: aA,T — 2 X olt and a\^^ — \pl X dnj where olt and are defined by equations 22 and 23, 
and aA — \/ ole X cka,!; X aA,T- 

^Tiie Peclet number is shown when provided by the author. In all cases the regions in the simulations for which parameter values 
are quoted were efficient convection, with AV < 10^^, and excluded the super-adiabatic layers in the surface convection models 
where parameters deviate significantly from those quoted here. 

*^The number of pressure scale heights spanned by the convcctivcly unstable region. 

"^The value in parentheses is for region spanning the entire computational domain, including the stable bounding layers with the 
other value referring to the convective region. 
^Model ob.Sd.B: additional details in Table 2. 

'^Chan & Sofia (1989): The range in a.T and olk,v is calculated according to their Table 1 for the nearly adiabatic portion of the 
simulation where 10"^ < AV < 10^^. 

^Kim et al. (1996); The coefficient olt is based on their Fig. 6. The coefficients aA,t) and oa are plotted in their Fig. 9 and the 
range quoted in the table above is for values at least one pressure scale height from the boundaries. 

^Robinson et al. (2004): only the correlation between radial velocity and temperature fluctuation is provided, which is a good 
surrogate for ctB- For the solar and subgiant cases see their Figs. 7-9. 

5(a)Porter & Woodward (2000): In this paper the values for a^, aT, and aA are quoted using the same definitions as in this 
study, (b) The lower value quoted by these authors for oa is a results of subtracting the kinetic energy flux from the enthalpy 
flux. The value in the parentheses is the mixing length oa according to the deflnition in note ^ above. 



TABLE 5 

Convective Boundary Layer Properties For Oxygen Shell Burning Models 



Model Time Interval rj 

(102 s) (10^ cm) 


h 

cm) 


(10* cm/s) 


(10* 


cm/s) 


(10^ 


]a 

cm/s) 


(10^ 


cm/s^) 


logE 


RtB 


ob.3d.B 


1.5, 2.7 


0.816 


1.287 


25.766±0.869 




0.6 




0.313 




0.574 


-1.095 


21.8 




2.7, 5.5 


0.842 


0.797 


8.252±0.180 








0.316 




0.966 


-1.616 


36.0 




5.5, 8.0 


0.861 


0.586 


5.171±0.179 








0.281 




1.062 


-1.789 


50.0 


ob.2d.c 


3.5, 5.7 


0.857 


0.191 


10.620±0.816 




0.9 




1.385 




1.422 


-2.154 


5.9 


ob.2d.C 


2.0, 4.0 


0.830 


1.776 


19.117±0.988 




0.5 




1.436 




1.010 


-1.887 


3.2 


ob.2d.e 


3.5, 8.0 


0.868 


1.900 


10.021±0.319 




2.5 




1.464 




1.334 


-2.289 


4.4 


ob.3d.B 
ob.2d.c 
ob.2d.C 
ob.2d.e 


3.0, 8.0 
3.5, 5.7 
2.0, 4.0 
3.5, 8.0 


0.429 
0.428 
0.430 
0.429 


0.057 
0.201 
0.193 
0.162 


-0.700±0.009 
-1.686±0.058 
-1.625±0.072 
-0.975±0.018 




0.50 
1.05 
0.65 
1.20 




0.479 
1.769 
1.434 
1.645 




30.686 
33.739 
32.160 
32.620 


-2.601 
-2.811 
-2.780 
-2.879 


418.6 
86.3 

101.7 
84.4 



^The rms fluctuations in the horizontal velocity at the interface location. 
'^The buoyancy jump across the interface. 



TABLE 6 

Convective Boundary Layer Properties For " Core Convection" Models 























Model Time Interval 




h 




Vexp 


(t[vh] 


Ab 


logB 


RtB 


(10^ s) 


(10" cm) 


(10^ cm) 


(10^ cm/s) 


(10^ cm/s) 


(10^ cm/s) 


(10^ cm/s^) 






msc.3d.B 


6.0, 10.0] 


1.374 


0.949 


1.754± 0.080 




2.011 


6.07 


-2.0594 


66 




10.0, 12.0] 


1.378 


0.897 


-0.020± 0.140 




1.878 


5.83 




72 




12.0, 15.0] 


1.382 


0.998 


2.731± 0.099 




2.411 


5.70 


-1.9459 


48 


insc.2(l.b 


6.0. 10.01 




1.319 


l.lOlit 0.390 




8.070 


6.13 


-2.7601 


9.2 



^The expansion velocity in these models remains very small with Vexp < 10 cm/s. 



